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ABSTRACT 

Context. The project Galactic Cold Cores has carried out Herschel photometric observations of interstellar clouds where the Planck 
satellite survey has located cold and compact clumps. The sources represent different stages of cloud evolution from starless clumps 
to protostellar cores and are located in different Galactic environments. 

Aims. We examine this sample of 116 Herschel fields to estimate the submillimetre dust opacity and to search for variations that 
might be attributed to the evolutionary stage of the sources or to environmental factors, including the location within the Galaxy. 
Methods. The submillimetre dust opacity was derived from Herschel data, and near-infrared observations of the reddening of back¬ 
ground stars are converted into near-infrared optical depth. We investigated the systematic errors affecting these parameters and used 
modelling to correct for the expected biases. The ratio of 250pm and J band opacities is correlated with the Galactic location and 
the star formation activity. We searched for local variations in the ratio r(250pm)/r(J) using the correlation plots and opacity ratio 
maps. 

Results. We find a median ratio ofr(250pm)/r(J) = (1.6 + 0.2) x 10 -3 , which is more than three times the mean value reported 
for the diffuse medium. Assuming an opacity spectral index J3 = 1.8 instead offi = 2.0, the value would be lower by ~ 30%. No 
significant systematic variation is detected with Galactocentric distance or with Galactic height. Examination of the r(250pm)/r(J) 
maps reveals six fields with clear indications of a local increase of submillimetre opacity of up to r(250/im)/r(/) ~ 4 X 10“ 3 towards 
the densest clumps. These are all nearby fields with spatially resolved clumps of high column density. 

Conclusions. We interpret the increase in the far-infrared opacity as a sign of grain growth in the densest and coldest regions of 
interstellar clouds. 

Key words. ISM: clouds - Infrared: ISM - Submillimeter: ISM - dust, extinction - Stars: formation - Stars: protostars 


1. Introduction 


The all-sky survey of the Planck satellite ( Tauber et al.|2010| has 
enabled a new approach to studying the earliest stages of star 
formation. The sub-millimetre measurements, with high sensi¬ 
tivity and an angular resolution down to ~ 4.5', have enabled 


* Planck (http://www.esa.int/Planck ) is a project of the European 
Space Agency - ESA - with instruments provided by two scientific 
consortia funded by ESA member states (in particular the lead coun¬ 
tries: France and Italy) with contributions from NASA (USA), and tele¬ 
scope reflectors provided in a collaboration between ESA and a scien¬ 
tific Consortium led and funded by Denmark. 

** Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im¬ 
portant participation from NASA. 

*** Tables 1 and E.l are only available in electronic form at the 
CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via 
http://cdsweb.u-strasbg.fr/cgi-bin/qcat? J/A+A/ 


detecting and classifying of a large number of cold and compact 
Galactic sources. These probably represent different phases in 
the evolution of dense interstellar clouds that leads to the for¬ 
mation of stars. Careful analysis of the Planck data has led to 
a list of more than 10000 objects that form the Cold Clump 


Catalogue of Planck Objects (C3PO, see Planck Collaboration 
|XXIII|2011 1 ) . At Planck resolution, it is not possible to resolve 
gravitationally bound cores even in the nearest molecular clouds. 
The low colour temperature of most of the sources (T < 14 K) 
strongly suggests that the Planck clumps must have high col¬ 
umn densities, possibly at scales not resolved by Planck , and 
they probably contain even dense cores. A significant fraction 
of the clumps may be transient structures produced by turbulent 
flows, however. 


Within the Herschel Open Time Key Programme Galactic 
Cold Cores , we have carried out dust continuum emission ob¬ 
servations of 116 fields that were selected based on Planck de- 
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tections listed in C3PO. The fields, which are typically ~40' 
in size, were mapped with Herschel PACS and SPIRE instru¬ 
ments ( [Pilbratt et al.||2010| Poglitsch et al.||2010[ |Griffin et~akj 
2010| ) at wavelen gths 100-500/ mi. The higher angular res olu- 
tion of Herschel ( Poglitsch et al.||2010 Griffin et al.||2Q 10 ) en¬ 
ables studying the internal structure of the Planck clumps, de¬ 
tecting individual cores, and, in conjunction with mid-infrared 
data, studying protostellar sources ( |Montillaud| [2014 sub mit¬ 
ted] ). The inclusion of far-infrared wavelengths helps to deter¬ 
mine the physical characteristics of the regions and, in particu¬ 
lar, to study the properties of dust emission. First results have 
been presented in |Planck Collaboration XXIII] (|2011|); |Planck| 
Coll aboration XXII| Q2011) and in |Juvela et al.| <20101 12011 


20121 (papers I, II, and III, respectively). Montillaud|( [20i4 sub^| 


mittedj ) presented the analysis of all clumps and cores found in 
our Herschel fields, including a comparison with the population 
of young stellar objects (YSOs). Further studies concentrated es¬ 
pecially on high latitude clouds ( [Malinen et al.| ( [Ml4| , Rivera- 
Ingraham et al. and Ristorcelli et al., in preparation). 

In this paper we concentrate on dust properties and es¬ 
pecially on the submillimetre dust opacity. Variations of dust 
emission properties have been investigated with far-infrared 
(FIR) and submillimetre observations of diffuse and molec¬ 
ular clouds, using data from IRAS, COBE, ISOPHOT, the 
PRONAOS balloon-borne experiment, and ground-based tele¬ 
scopes. It was shown that the low temperatures found in a sample 
of m olecular clouds ( [Laureijs et al. || 1 99T] |Abergel et al.||l 994 


et al 


1996) and the translucent Polaris Flare cirrus cloud (Bernarc 
[1999| cannot be explained by the extinction of the radi 


ation field. An increase of the dust emissivity by a factor of 3 
compared to the standard diffuse value was needed to reproduce 
the cold temperatures observed in the Taurus filament LI506 


shown an opacity increase by a factor between 2 to 4 ( 

Cambresy 

et al. 

2001 

Kramer et al. 2003; Bianchi et al. 2003j; 

del Burgo 

et al. 

2003 

Kiss et al. 2006; Ridderstad et al. 2006 

Lehtinen 

et al. 

2007 

). More recently, similar results have been obtained 


Herschel observations of the L1506 filament, Ysa rd et aL] ( |2013| ) 
characterised the dust evolution toward the dense part of the fil¬ 
ament. The dust emissivity in the outer layers of the filament 
was found to be consistent with standard grains from the dif¬ 
fuse medium, whereas the emissivity increases by a factor of ~2 
above gas densities of a few times 10 3 cm -3 . This change has 
been attributed to the formation of fluffy aggregates in the dense 
medium, resulting from grain coagulation, as suggested by pre¬ 


vious studies ([Cambresy et al. 2001 

Stepnik et al. 2003; Bernard 

et al. 1999 

del Burgo et al. 2003; 

Kiss et al. 2006; Ridderstad 

et al. 2006 

). The average size of aggregates required to fit the 


shape, structure, material composition, and accretion of mantles 
can also all contribute. 

Moreover , Malinen et al.| ( |2011| ); | Juvela & Ysard| ( |2Q12| ); 
Ysard et al. ( |2012| ) have investigated the impact of radiative 


([Juvela et al. |2011[ |Planck Collaboration XXV||2011[ |Martin| 
et al. 2012| Roy et al. 2013 ). In their det ailed modelling of 


transfer on the results derived from observations under the as¬ 
sumption of a single average colour temperature. They showed 
that the mixing of different temperatures along the line of sight 
produces a tendency that is opposite to the observed one. They 
concluded that in the absence of internal heating sources, the 
observed emissivity increase toward dense clouds cannot be ex¬ 
plained by radiative transfer effects. It must originate in intrinsic 
variations of the optical properties of the grains. 

It is, however, important to note that the dust emissivity in¬ 
crease is not systematically observed in the inters t ellar medium 
(ISM, see jNutter et al.| ( |2()08| ); | Juvela et al.| ( [20Q9| ); |Paradis et alT] 
( 2009])) . These intriguing results call for a broader investiga¬ 
tion, making use of the large observations statistics provided by 
Herschel and Planck , probing different Galactic environments. 
The key questions are still open today: when, where, and how 
dust evolves between diffuse and dense regions, what the phys¬ 
ical conditions enhancing (or preventing) the efficiency of the 
coagulation process are, what the time scales are, and whether 
the process is directly related to specific stages in the cloud or 
core evolution. Understanding these questions is critical since 
knowing the dust opacity has a direct impact on many key pa¬ 
rameters derived from dust emission, such as the column densi¬ 
ties, masses, and volume densities of the clouds. For this reason, 
it is also necessary to investigate the possible systematic effects 
on the emission and extinction measurements that could cause 
errors in the opacity estimates. 

The structure of the paper is as follows: The observations are 
described in Sect. [2] The main results are presented in Sect. [3] 
including the estimates of submillimetre and near-infrared (NIR) 
optical depths, the correlations between these variables, and 
the correlations with environmental factors. The results are dis¬ 
cussed in Sect.[4]before we list the final conclusions in Sect. [5] 


2. Observations and basic data analysis 

2.1. Target selection 

The selection of the Herschel fields is described in Juvela et al. 


(2012) and an overview of all the maps is given in Montillaud 


FIR, submillimetre, and extinction profiles in the LI506 filament 
is about 0.4 /i m ( [Ysard et al.||2013| ). This value is close to the 
smallest grain size needed to scatter light efficiently in the mid- 
IR, which produces the ’coreshine’ observed toward a number of 
dense cores, which has also been interpreted as a result of grain 
growth ( [Pagani et al.|2010t|Steinacker et al.|2010] ). 


On the theoretical side, various dust optical property calcu¬ 
lations have predicted a significant increase in the emissivity of 
aggregat es at long wavelengths compared to compact spherical 
grains ([Wright |1987[ |Bazell & Dwek||1990[ |Ossenkopf| 1993[ 


Ossenkopf & Henning 1994; Stognienko et al. 1995; Kohler 
et al.|201 1||2012| ). This variation is shown to be mainly due to the 
increase in the porosity fraction with aggregate growth, but the 


(2014 submitted). We only repeat the main points here. 

Planck sub-millimetre observations, together with IRAS 
lOOyurn data, enabled the detection of over 10000 compact 
sources in which the dust is significantly colder than in the sur¬ 
rounding regions ( [Planck Collaboration XXIH]|2011[ ). The de¬ 
tection procedure is based on this colour temperature difference 
and, furthermore, limits_the size of the detected clumps to values 
below ~12' ( Mon tier et al.|2 010). The sources are believed to be 
Galactic cold clumps or, at larger distance, entire clouds ([Planck 
[Collaboration XXIII|2011| ). 

The fields for Herschel follow-up observations were se¬ 
lected using a binning of Planck cold clumps with respect to 
the Galactic longitude and latitude, the estimated dust colour 
temperature, and the clump mass. At the time of source selec¬ 
tion, distance estimates existed for approximately one third of 
the sources in C3PO and, therefore, some sources of unknown 
mass were also included. The binning ensured full coverage 
of the clump parameter space, especially of the high Galactic 

o 

latitudes and of the outer Galaxy. Galactic latitudes \b\ < 1 
were excluded because that area is covered by the Hi-GAL pro¬ 
gramme ( Molinari et al.|[2010[ ). Similarly, regions included in 
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othe r Herschel key pr ogrammes such as the Gould Bel t sur¬ 
vey ( [Andre et al.||2010) and HOBYS ( Motte et al.||2010| ) were 
avoided. 

A total of 116 separate fields were observed. The SPIRE 
maps are on average over 40' in linear size, with an average 
area of ~ 1800 (arcmin) 2 . The PACS maps are smaller, with 
an average area of ~ 660 (arcmin) 2 . Most fields contain more 
than one Planck clump, the maps altogether covering ~350 indi¬ 
vidual Planck detections. The range of probed column densities 
extends from diffuse fields with /V(H 2 ) ~ 10 21 cm -2 to cores 
with A(H 2 ) > 10 23 cm -2 . The fields are listed in Table [1] and the 
Herschel observation numbers are included in Table! 


polated Planck and IRIS data. The zero points were calculated 
iteratively, including colour corrections calculated using colour 
temperatures that were estimated from SPIRE data with a fixed 
spectral index of J3=2.0. 


2.2.3. Calculating submillimetre optical depth 

The Herschel maps were converted into estimates of dust optical 
depth at 250 jam. The surface brightness maps were convolved to 
a common resolution of 40" , and colour temperatures were cal¬ 
culated by fitting the spectral energy distributions (SEDs) with 
modified blackbody curves with a constant opacity spectral in¬ 
dex of p =2.0. The 250/im optical depth was obtained from 


2.2. Herschel data 

2.2.1 . Herschel data reduction 


r(250//m) = 


7 y (250/rm) 

By(T) 


( 1 ) 


The fields were mapped with the SPIRE instrument at wave¬ 
lengths 250, 350, and 500 jam and with the PACS instrument 
at wavelengths 100 and 160/im. One field was observed with 
SPIRE alone (G206.33-25.94, part of the Witch Head Nebula, 
IC 2118). The Hersch el ob servations are discussed in detail 
in ( [Juvela et al.||2012| ) and ( |Montillaud||2014 submitted] ). The 
SPIRE observations at 250/mi, 350/rni, and 500 jam were re¬ 
duced with the Herschel Interactive Processing Environment 
HIPE v. 10.0.0, using the official pipeline with the iterative 
destriper and the extended emission calibration options. The 
maps were produced with the naive map-making routine. The 
PACS data at lOO/tm and 160yum were processed with HIPE 
v. 10.0.0 up to level 1, and the maps were then produced with 
Scanamorphos v20 (Roussel |2013| . In the order of increasing 
wavelength, the resolution of the maps is 7", 12", 18", 25", 
and 36" for the five bands. The raw and pipeline-reduced data 
are available via the Herschel Science Archive, the user-reduced 
maps are available via ESA sit^] 

The accuracy of the absolute calibration of the SPIRE ob¬ 
servations is expected to be better than For PACS we as¬ 
sume a calibration uncertainty of 15%. This is a conservative 
estimate that is compatible with the differences of PACS and 
Spitzer MIPS measurements of extended emission^] 

2.2.2. Estimating intensity zero points 

To determine the zero point of the intensity scale, we compared 
the Herschel maps with Planck data complemented with the 
IRIS version of the IRAS lOO/rni data ( [Miville-Deschenes ~&| 
|Lagache||2005| ). The Planck and IRIS measurements were in¬ 
terpolated to Herschel wavelengths using fitted modified black- 
body curves, B v (T du st )v^, with a fixed value of the spectral in¬ 
dex, p = 2.0. The linear correlations between Herschel and the 
reference data were extrapolated to zero Planck (+IRIS) sur¬ 
face brightness to determine offsets for the Herschel maps. For 
SPIRE the uncertainties of these fits are typically -lMJysr -1 
at 250 yum and at longer wavelengths smaller in absolute value. 
The derived intensity zero points are independent of Planck cal¬ 
ibration and of any multiplicative errors in the comparison. For 
PACS the correlations are often less well defined, and the zero 
points were set directly based on the comparison of the aver¬ 
age values of the Herschel maps and the corresponding inter¬ 


using the fitted 250 yum intensity 7 y (250/im) and the colour tem¬ 
perature T . The calculations were made with 250-500 /am data 
and 160-500//m data. The fits were weighted according to 15% 
and 7% error estimates for PACS and SPIRE surface brightness 
measurements, respectively (see Sect. |2. 2. 1] ). 

The assumed value of p= 2.0 may be appropriate for dense 
clumps, although at lower column densities the average value 
is lowe r, p ~ 1.8 ( [Boulanger et al. 1996 [ [Planck Collaboration 


XXV 


201 1| ), and the value of p may further depend on 


the G alactic l ocation and the wavelength range (e.g. |Planck| 
[Collaboration Int. XIV|[2014) . If the true value of p were 1.8 
instead of 2.0, our colour temperature estimates would be higher 
by ~1 K and the r(250yum) values lower by ~30%. Furthermore, 
if the values of p were correlated with column density, the slope 
of r(250yum) vs. T/ would be similarly affected. We return to 
these effects in Sects. 13.7l andl4l 

To estimate the statistical uncertainty of r(250yum) values, 
we used Markov chain Monte Carlo (MCMC) runs. The prior 
distribution of temperature values is flat but limited between 
5.0 K and 35 K. In addition to the relative errors quoted above, 
we included the uncertainty of the intensity zero points. These 
are typically much smaller than the assumed relative errors, but 
may be important at low column densities, especially at 160 p m. 
The zero-point errors are systematic but are included simply as 
another component of statistical noise. Their effect is thus re¬ 
flected in the error estimates of individual pixels. The error dis¬ 
tribution of r(250/im) is nearly Gaussian, and we used the stan¬ 
dard deviation of the MCMC r(250/mi) samples as the error esti¬ 
mates. These estimates were calculated separately for each pixel 
of the r(250yum) maps. 

Because of line-of-sight temperature variations, the derived 
r(250yum) e stimates probably systematically underestimate the 
true values ( [Shetty et al.|2009[|Malinen et al.|201l| ). We cannot 
directly determine the magnitude of these errors but, with some 
assumptions, we can use radiative transfer modelling to estimate 
the magnitude of the bias. The simulations, described in detail 
in Appendix [Cj were used to derive bias maps that are taken into 
account when the data were correlated with tj values. 

2.3. Near-infrared data 


We used the Two Micron All Sky Survey (2MASS, [Skrutskie 


1 http://herschel. esac. esa. int/UserReducedData. shtml 

2 SPIRE Observer’s manual, 

http://herschel. esac. esa. int/Documentation. shtml 
5 http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb 


et al.||2006| ) to derive estimates of dust column density that 


are independent of dust emission. We used the method NICER 
(Lombardi & Alves 2001) and the standard extinction curve 
( Cardelli etal. 19 89] ) to convert the reddening of the background 
stars to estimates of J-band optical depth, tj. Because the cal- 
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culations involve only near-infrared bands, the results are ex¬ 
pected to be insensitive to the value of the ratio of total to se¬ 
lective extinction, Ry ( [Cardelli et al.||1989| . The shape of the 
NIR extinction curve is believ ed to be relatively stable, even 
at high extinctions (e .g. Draine 2003a; Indebetouw et al. 2 005 [ 


Roman-Zuniga et al.|2007[ Ascenso et al.| 


Lombardi et al. 


2013[|Wang & Jiang|2014| ). Some variations are observed with 

Galactic location and/or density, but generally only at a level of 
5% of the N IR power-law index (e.g. |Stead & Hoare 2009[|Fntz| 


et al. 2011). This question is discussed in more detail in Sect. 4.2 


The t(J) values are derived using both the J-H and H-K colours 
but, with the extinction curve used, we have the correspondence 
of E j- K = 0.65 r(J). Flags in the 2MASS catalogue were used 
to avoid galaxies (extJkey not null or gaLcontam not zero) and 
sources with uncertain photometry (ph.qual worse than C). 

Five of our fields are fully covered in the VISTA Hemisphere 
Survey, VHS (McMahon et al. |2013| ), which has more than ten 
times the sensitivity of 2MASS (in H band VHS has a 5-cr de¬ 
tection threshold of 19.0, compared to a 2MASS point source 
catalogue completeness limit of ~16 mag). One of these fields 
is too distant to obtain a reliable extinction map, but the data for 
the four other fields were analysed and the results compared with 
those obtained with 2MASS data. The fields are G4.18+35.79 
(LDN 134), G21.26+12.11, G24.40+4.68, and G358.96+36.75. 
For the fields G21.26+12.11 and G24.40+4.68, only J- and 
Ks-band data exist. The data are available in VISTA Survey 
Archiv^] , and the VISTA Data Flow System pipeline process¬ 
ing and science archive are described in |Irwin et al.| ( 2004) and 
|Hambly et al.| ( |2008 ). 

Extinction maps are produced by averaging extinction es¬ 
timates of individual stars with a Gaussian weighting func¬ 
tion with FWHM=180". We also tested a higher resolution of 
FWHM=120". For distant sources, the extinction of the target 
clouds cannot be reliably reproduced because of the poor res¬ 
olution and the increasing number of foreground stars. This is 
the main factor that limits the number of fields where the ratio 
r(250/im)/r(/) can be reliably estimated. The extinction mea¬ 
surements can be significantly biased even in nearby fields if 
these contain steep column density variations. No special steps 
were take n to eliminate the con tamination by foreground stars 
(see, e.g., Schneider et al.|20n) , apart from the sigma clipping 
procedure that is part of the NICER method and was performed 
at 3 cr level. The reliability of the extinction maps and the bias 
caused by sampling problems and the presence of foreground 
stars was examined with simulations (see Sect. 0 . The results 
of these simulations are used to derive maps of the expected un¬ 
certainty and the bias of the r(/) values for each field. 


2.4. Correlations between sub-millimetre and NIR opacity 

The ratio k of sub-millimetre opacity r(250/im) and the NIR 
opacity tj was estimated for all 116 fields. The r(250yum) 
maps were convolved to the 3' resolution of the tj maps. The 
r(250yum) and the r(J) data were read at 90" steps (half-beam 
sampling), excluding the map borders where the result of the 
convolution to 3' resolution is poorly defined. For local back¬ 
ground subtraction, only areas where the signal was more than 
2 cr above the average value of the reference area were used 
(see |2.2.2| ). Here cr is the standard deviation of the values in the 
reference region. This is a conservative limit because part of the 
fluctuations is caused by real surface brightness variations and 
not by noise alone. 

4 http://homs.roe.ac.uk/vsa/index.html 


For r j the error estimates were provided by the NICER rou¬ 
tine. For r(250//m) these were obtained from MCMC calcula¬ 
tions (see Sect. |2.2.3| ). The comparison between the different 
cases (for example, regarding the use of 160/im data, back¬ 
ground subtraction, or gradient corrections) provides informa¬ 
tion on the uncertainty caused by some sources of systematic 
errors. 

The r(250yum) vs. r(/) points of individual fields and sam¬ 
ples of fields were fitted with a linear model to derive the ratio 
r(250/mi) /r(J). These total least-squares fits take into account 
the uncertainties in both variables. The fits were made using ei¬ 
ther all data points or only data below or above a given r(/) 
limit. To reduce the bias caused by these cuts, the data were di¬ 
vided with the help of a preliminary linear fit to all data points 
(see Sect. |3.1| for details). The limiting value of r(/) thus cor¬ 
responds to a position on this line, and the cut was performed 
using a line that is perpendicular in a coordinate system where 
the average uncertainties of the two variables are equal. 

The r(250/im)/r(/) ratios were also calculated for alterna¬ 
tive versions of the r(250yum) data, using local background sub¬ 
traction or using ancillary data in an attempt to correct for possi¬ 
ble large-scale errors in the surface brightness data. These alter¬ 
native data are discussed in Appendix \K\ 


3. Results 

3.1. Apparent r(250pm)/r(J) values 

We calculated t(J) and r(25 0yum) maps of the 116 fields as de¬ 
scribed in Sects. 1^2] and |2. 3 1 The correlations between r(J) and 
r(250yum) were calculated at a resolution of 180". In addition 
to the full range of column densities, the relationships were ex- 
amin ed s eparately below and above the limit of r(/)=0.6 (see 
Sect. |2.4| ). This corresponds to visual extinctions Ay ~ 2.3 mag 
and Ay ~ 2.0 mag for the Ry values of 3.1 and 5.0, respectively 
( [Cardelli et al.|1989] ). Instead of a higher limit, we selected the 
relatively low number of r(/) = 0.6 to maximise the number of 
fields where a linear fit could also be made above the r(/) thresh¬ 
old. The r(250yum) values were derived from Herschel data with 
either 250-500 pm or 160-500 pm (see Appendix |A| for analysis 
with additional alternative data sets). 

In a given field, the number of points either below or above 
the t(J) limit is often insufficient to determine any reliable value 
for the slope k = Ar(250yum)/Ar(/). In a few fields no reliable 
value of k can be determined at all, mainly because of the low 
quality of the r(J) data. This especially affects the most dis¬ 
tant fields because of the contamination by foreground stars and 
because the structures are too small to be resolved with the 3' 
beam. The formal errors of the k parameter were used to exclude 
the clearly unreliable fits. The criterion 6k/k <0.1 leaves in the 
default case 106 fits to all data in a field, 103 fits below r(/)=0.6, 
and 38 fits above r(/)=0.6. These fits appear relatively reliable 
also based on visual inspection. 

Figure [T] shows an example of the recovered dependence be¬ 
tween t(J) and r(250yum) values, including linear fits to the three 
r(/) ranges. In this example, the slope appears to become steeper 
as t(J) increases. This might be an indication of an increase in 
the dust submillimetre opacity, which in turn might be attributed 
to grain growth (e.gJOssenkopf & Henning 1994; Stepnik et al. 
|2003[ |Ormel et~aL||2011HYsard et al.||2013| ). However, before 
drawing any such conclusions, we must consider the systematic 
effects that affect the two parameters. Figure [2] shows a summary 
of all the Ar(250pm)/Ar(J) values where r(250/im) values are 
based on Herschel 250-500 pm data. Before any bias corrections 
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Fig. 1. Relation between r(250yum) and r(/) in the field 
G95.76+8.17. The black solid line is a linear weighted total 
least-squares fit to all data points. The blue and red points and 
lines of the corresponding colour show the data and the fits be¬ 
low and above the threshold of r(/) =0.6, the dashed line in¬ 
dicates the division. The values of the slopes k are given in the 
plot. Error bars are shown for a set of random data points. 


(see below), the values are seen to cluster around ~ 2.0 x 10 -3 , 
with some tendency for higher values in the higher r(/) range. 

Figure [3] summarises the statistics of the dust opacity mea¬ 
surements as histograms, including all fits where the formal error 
of the slope of the least-squares fit r(250yum) vs. r(/) is below 
10 %. 


We need to observe a sufficient number of background stars 
for each resolution element, even at high column densities. This 
means that r(/) estimates and the comparison with submil¬ 
limetre emission can be made only at a low resolution (2-3'). 
Averaged over such large areas, th e sta tistical uncertainty of 
Herschel data is very small. In Sect. |3.2| we show that the bias 
is probably also dominated by the errors in r(/). In the follow¬ 
ing, we rely mainly on the Herschel data set that consists of ob¬ 
servations 250-500yum (the “default” data set). There are three 
reasons. First, in theory, the inclusion of the 160yum data re¬ 
duces statistical uncertainty of the colour temperature estimates 
but increases systematic errors caused by line-of-sight tempera¬ 
ture variations ( [Shetty et al.|2009|[Malinen et al.|201 1| ). Second, 
because of the smaller (and, for parallel mode, different) area 
covered by the PACS observations, the use of the 160yum band 
significantly reduces the area where correlations with r(/) can 
be calculated. Third, 160yum data may be affected by additional 
systematic effects related to the relative calibration of the two 
instruments, uncertainties in the zero-point determination (inter¬ 
polation between IRAS and Planck channels and the contribu¬ 
tion of stochastically heated grains in the IRAS 100/mi band) 
and to imperfections in the map making that coul d be i ncreased 
by the smaller size of the PACS maps (see Sect. A.2 ). We are 
particularly interested in the coldest regions where 250-500 /i m 
data provide adequate constraints on the dust temperature. 


3.2. Bias in t(J) values 

Bias in r(J ) values is very likely a significant problem, espe¬ 
cially for distant fields in which all high column density struc¬ 
tures are not resolved and the results begin to be affected by 
foreground stars. Both effects decrease the r(J) estimates, es¬ 



Fig.3. Comparison of Ar(250yum)/Ar(/) values in three r(J) 
intervals (three frames), obtained using either three of four 
Herschel bands in deriving the r(250yum) values. The two his¬ 
tograms without hatching (thick outlines) show all fields where 
the estimated uncertainty is below 10%. The two hatched his¬ 
tograms contain only the intersection with better than 10% ac¬ 
curacy with both three and four bands (79, 83, and 14 fields for 
the three panels, respectively). 


pecially towards column density peaks. We estimated the ex¬ 
tent of the problem with simulations using the stellar statistics 
in low-extinction areas near each field. The contamination by 
foreground stars was evaluated with the help of the Besangon 
model of the Galactic stellar distribution ( [Robin et al.|2003| ). We 
used the Herschel column density maps as a model of the col¬ 
umn density structure, simulated the distribution of foreground 
and background stars, analysed the simulated observations with 
NICER routine, and compared the results with the known input 
t(J) map. The procedure is described in detail in Appendix [B] 
We obtained for each field a map of the expected systematic rel¬ 
ative error in t(J) that gives a multiplicative correction factor 
for the r(/). The simulations do not consider the effect of cloud 
structures at scales below 18" but the procedure probably pro¬ 
vides a reasonable estimate of the magnitude of the effect. 

We repeated the analysis of the previous section and replaced 
the original r(/) maps with bias-corrected estimates. Figure [4] 
compares the Ar(250/im)/AT(/) distributions for the default case 
with and without bias correction. The statistics include all fits 
for which the formal error of the slope is below 10%. This cor¬ 
responds to Fig. [3] but the number of points is different. Because 
the bias corrections depend on the cloud distance, fields without 
distance estimates had to be dropped. However, in the r(/) > 0.6 
interval the number of fields fulfilling the 6k/k <0.1 criterion 
has doubled. 
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Fig. 2. Slopes k = Ar(250yum)/Ar(/) for all cases with uncertainties 5k/k <0.1. The black, blue, and red symbols correspond to 
values derived for the full r(/) range and for data below and above the limit of r(/)=0.6. The values of r(250yum) have been derived 
from SPIRE data without the subtraction of the local background. Neither r(250yum) nor r(/) has been corrected for the expected 
bias. 



Fig. 4. Comparison of Ar(250/dm)/Ar(J) distributions without 
bias corrections (“Default”) and with bias corrections applied 
either to r(J) or to both r(J) and r(250/im). The three frames 
correspond to different ranges of r(J) values. 


3.3. Bias inr( 250/rni) values 

The systematic errors in r(250yum) values were estimated with 
radiative transfer modelling. The line-of-sight temperature vari¬ 
ations are expected to be the main source of error that, in stan¬ 
dard analysis, leads to overestimation of the mass-averaged dust 


temperature and, subsequently, to underestimation of r(250yum) 
(e.g. Ysard et al.|2012) . 


We constructed for each field a three-dimensional radia¬ 
tive transfer model that covered a projected area of 30' x 30' 
with a 10" pixel size. The modelling assumed spatially constant 
dust properties, the dust model ( |Draine|2003b| ) corresponding to 
Rv=5.5 (see Appendix [C] for details). The density distribution 
and external heating were adjusted until the model exactly re¬ 
produced the observed 350/im surface brightness and, for the 
area above median column density, the average 250yum/500 jum 
ratio. The model-predicted surface brightness maps were anal¬ 
ysed as in the case of the actual observations, to produce maps 
of r(250yum). To estimate the bias, these values were compared 
to the actual r(250yum) values of the model to derive multiplica¬ 
tive correction factors. 


The results depend on th e ass umed cloud structure in the 
line-of-sight direction (Juvela et alJ 2013|. In our models, the 


line-of-sight density distribution only has one peak. This en¬ 
hances temperature contrasts and increases our bias estimates. 
On the other hand, the densest observed cores are probably even 
more compact, and in their case we may be underestimating the 
bias. If the clouds contain embedded sources, the actual bias may 
again locally be very different and often lower than predicted by 
our models. Although the bias estimation is more difficult than 
for r(/), the models should again provide a reasonable estimate 
of the magnitude of the effect. The relative systematic errors are 
smaller in r(250yum) than in r(J) so that their effect on the final 
result is less strong. 

The k = Ar(250yum)/Ar(/) values were re-calculated includ¬ 
ing bias corrections in both variables. The resulting histograms 
are included in Fig. [4] The bias correction makes the distribu¬ 
tions significantly narrower. Figure [4] also shows that the correc¬ 
tions are much stronger for r(J) than for r(250/im). As a result, 
the average value of Ar(250yum)/Ar(/) now decreases with in¬ 
creasing t(J). The number of fields fulfilling the 5k/k <0.1 cri¬ 
terion has doubled to 76 fields (using SPIRE bands). Compared 
to the original data, the median values of k x 10 4 have decreased 
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Fig. 5. Slope values k = r(250/im)/r(/) of Fig. [ 2 ] as a function of estimated distance (frames a, b ), galactocentric distance (frames 
c, d ), and galactic height (frames e , /). The left frames show the original slopes, the right frames the slopes after the bias corrections 
of r(250yum) and r(J). The black, blue, and red colours correspond to the full r(J) range and to data below and above r(/)=0.6, 
respectively. The solid curves with the same colours are the weighted moving averages (window sizes 30% in distance, 800 pc in 
Galactocentric distance, and 100 pc in Galactic height). All r(250/mi) values are calculated with SPIRE bands alone. 


from 20.2, 21.1, and 23.4 to 15.3, 16.0, and 12.2, the numbers 
corresponding to the full r(J) range, data below r(J) = 0.6, and 
data above r(J) = 0.6, respectively. The strong change in the k 
values suggests that the uncertainty of k is probably often several 
tens of per cent, especially in the r(J) > 0.6 interval. Therefore, 
Fig.[4]does not exclude a systematic increase of k as the function 
of r(/), if that becomes visible only at high column densities. 

Figure [5] displays the slopes k = Ar(250yum)/Ar(/) as func¬ 
tion of distance and Galactic location. The left frames show the 
relations without and the right frames with the bias corrections 
applied to r(J) and r(250yum). Only fits with 5k/k <0.1 are in¬ 
cluded. The original data showed some trends, including an in¬ 
crease in k as a function of distance and galactocentric distance. 
The first is visible especially in the high r(J) interval, but was 
expected because r(J) values of distant sources can be severely 
underestimated. After bias corrections, the scatter of the k values 
is significantly reduced. This suggests that the corrections are of 
correct magnitude. The distance dependence has changed so that 


in the corrected data there is a slight decrease in k values as a 
function of distance. This could point to some over-correction 
of the t(J) estimates, although the bias correction should not 
only depend on distance, but even mainly on the cloud structure. 
However, the decrease of k can be an indication of selection ef¬ 
fects or direct resolution effects. For example, higher k values 
might be found in individual dense clumps that are only resolved 
at short distances. 

There is little difference between the k values found in the 
three t(J) intervals. In the next section we examine in more de¬ 
tail the global t(J) dependence of the r(250yum)/r(/) ratios, es¬ 
pecially regarding the highest observed column densities. 

3.4. Global relation r(250/mi) vs. r(J) 

To further test the hypothesis that r(250yum)/r(/) ratios may 
change systematically as a function of column density, we car¬ 
ried out non-linear fits r(250/dm) = A + Bx r(J) + Cx r(/) 2 . 
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The fits were first performed using the combined data of all 
fields. To reduce the mismatch in the zero levels of individual 
fields, we subtracted from each r(/) and r(250yum) map the lo¬ 
cal background using the off regions listed in Table [I] The off re¬ 
gions are not completely void of emission but provide a common 
reference point for the quantities. Thus, the relation is expected 
to develop via the origin for each field separately, the parameter 
A being close to zero for the combined data as well. The sign of 
the fitted parameter C indicates the possible increase or decrease 
of k- Ar(250/im)/Ar(/) as the function of column density. 

Figure [6] shows the results obtained with the bias-corrected 
data. We included data from all fields in which individual lin¬ 
ear fits had Sk/k < 0.2, thus relaxing the previous constraint 
of Sk/k < 0.1. The second-order polynomial was fitted to all 
data and separately to data points with r(/) > 1.0. In the previ¬ 
ous section a threshold value of r(/) = 0.6 was used. However, 
some 70% of all points are below r(/) = 0.6 and, when included, 
they dominate the fits that systematically underestimate the data 
above r(/) ~ 5. With the combined data set, there are enough 
high column density data points so that the lower limit can be 
moved upwards. The use of the r(/) = 1.0 threshold enables an 
adequate fit to all data with higher r(/) values. The r(/) calcu¬ 
lations employ a different off region for each field. These may 
contain different amounts of extinction, which leads to small 
relative shifts along the r(/) axis. Based on dust emission, the 
extinction in the off regions is typically r(/)=0.2-0.4. The un¬ 
certainty of the relative zero points contributes to the scatter in 
Fig. [6] but the effect is weaker than the total dispersion and the 
non-linearity seen at high extinctions. 

To prevent the r(J) cut itself from biasing the fits, the data 
were selected using lines perpendicular to a linear least-squares 
line fitted to all data (cf. Sect. |2.4| ). Thus, the quoted r(/) limits 
correspond to a point on the fitted line, and the cut itself is per¬ 
pendicular to the fitted line. All fits take into account the uncer¬ 
tainties in both variables, which are assumed to be uncorrelated. 
The error distributions of the parameters A-C were calculated 
with an MCMC. 

The sign of the parameter C depends on the range of r(/) 
values but is less dependent on the field selection, for example 
regarding the Sk/k limit that was used to select the fields. Most 
fields with high r(/) values (and thus with a wide dynamical 
range) have low values of Sk/k. When all points are included, 
the value of the parameter C is negative, but the fit is very poor 
at high r(/). When pixels r(/) <1.0 are excluded, C becomes 
positive (see Fig. [6^), which points to an increase in the sub¬ 
millimetre opacity above r(/) ~ 1. Systematic additive errors 
in either parameter might also explain the different behaviour at 
very low r(/). When the fit is made using all data r(/) > 0.6 
(not shown), the parameter C is marginally positive, but beyond 
r(/) = 10 the fitted line is below all the data points. The second 
frame of Fig. [6] shows the fits when pixels with colour tempera¬ 
tures above 14 K are excluded. The values of C are now higher 
and positive even when data r(/) <1.0 are not excluded. The 
best fit to the high r(J) end of the relation is still obtained by 
excluding data with r(/) < 1.0, this results in the relation 

r(250jum) = 0.73 xl0" 3 +1.25 xlO" 3 r(/) + 0.11 x 10“ 3 t(/) 2 .(2) 

The formal error estimates of the parameters A-C are of the order 
of 5%, probably lower than the systematic uncertainties. All data 
beyond r(/) ~ 5 are affected by large bias corrections and, con¬ 
sequently, the value of C also depends on the accuracy of these 
corrections. Thus, Fig. [6] strongly suggests but does not yet pro¬ 
vide a final proof of the variations of the ratio r(250/im)/r(/) . 
The positive offset A=0.7xl0 -3 results from the facts that at low 



Fig. 6. Fit of r(250//m) = A + B x r(J) + C x r(/) 2 to the com¬ 
bined data of all fields in which individual linear fits showed a 
strong correlation with Sk/k < 0.2. The blue and the green lines 
correspond to fits to the full column density range and to data 
points t(J) >1.0 alone, respectively. In the second frame, only 
data with colour temperatures below 14 K are used. 
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Fig. 7. Distributions of the parameters of the fit r( 250/rm) = 
A + B x r(/) + C x r(/) 2 . The fit is limited to data with r(/) > 1. 


column densities the relation is linear, the curvature increases 
only beyond r(/) ~ 5, and the lowest data points r(/) <1.0 are 
not part of the fit. 


Figure[7]shows the error distributions of the parameters A-C. 
The fit was made to data r(J) >1.0 for all fields with a linear fit 
accuracy Sk/k < 0.2. In most fields, the formal error estimates 
of St(J) and Sr( 250/rni) are smaller than the actual scatter of 
points. Therefore we used the residuals of the linear fits before 
the MCMC calculation to determine a scaling factor, typically 
2.0-3.0, that makes the error estimates in each field consistent 
with the actual scatter. Even after this increase of uncertainties, 
MCMC gives a 100% probability for a positive value of C. In 
reality, the result is not that strong because the uncertainty may 
be dominated by systematic errors. The sign of C was already 
seen to change depending on the range of r(/) values fitted. The 
result also depends on a relatively small number of fields with 
data above r(J) > 5. Therefore, we must consider the r(250/rm) 
vs. r(/) relation in individual fields in more detail. 
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3.5. Correlations in selected fields 


Figure [6] showed hints of an increase of the r(250yum)/r(/) val¬ 
ues as the column density increases, but the global statistics may 
be confused by the mix of different fields. Furthermore, the sign 
of the C parameter is determined by the highest r(/) points that 
originate in a small number of individual fields. Each field may 
be affected by different systematic effects related to the surface 
brightness zero points, distance uncertainty (via bias correction), 
and differences in the local radiation field. Several diffuse fields 
are even entirely below r(J) ~ 1.0. Therefore, we also need to 
examine the fields individually. 

Three criteria were used to select a subset of fields. We re¬ 
quired that (1) the uncertainty of the fitted parameter C is below 
0.3 x 10 -4 , (2) there are at least ten data points (selected from 
the maps at 90" steps) with r(J) above 0.6, and (3) the bias cor¬ 
rections change the slope of the linear fit of r(250/mi) vs. r(J) by 
less than 30%. The first two criteria ensure that there are enough 
data points at large r(J) with a small scatter to gain some insight 
about the column density dependence of the Ar(250yum)/Ar(/). 
The third criterion excludes distant fields for which the uncer¬ 
tainty of the bias correction of r(J) renders the results uncertain, 
even for the apparently well-defined relation between r(/) and 
r(250yum). The selection leaves 23 fields with distances mostly 
in the range 100-500 pc. There are two exceptions, G216.76-2.58 
and G111.41-2.95, for which the estimated distances are 2.4 kpc 
and 3.0 kpc. We kept the two fields in the sample even though the 
results are known to be unreliable because of the large distance. 

To avoid underestimating the fit errors, we scaled the error 
estimates of r(J) and r(250/rm) up to correspond to the actual 
scatter of points (see Sect. 3.4). All observations, sampled at 90" 
steps, are fitted with a linear model r(250/im) = b + k x r(J) us¬ 
ing total least-squares. We continued to use as the default data 
set one with r(250yum) derived from SPIRE bands alone, includ¬ 
ing bias corrections in r(250/im) and r(J). However, for com¬ 
parison, we also examined results obtained with four Herschel 
bands (160-500 yum; including the bias corrections) and, finally, 
with three Herschel bands but without any bias corrections. 

The non-linear fits were made using MCMC (with 2 x 10 5 
samples per field) and bootstrap sampling (2000 realisations per 
field). Both fits used total least-squareq^Jand the error estimates 
of the individual points. The parameter values and the uncertain¬ 
ties derived with these two methods should be similar, except 
for rare cases in which the result depends on a small number of 
influential points, which are always present in the MCMC cal¬ 
culation, but not in all bootstrap samples. 

Figures [8p| show the results of these fits. Each frame shows 
the values of the linear slopes for the three cases discussed 
above. The parameters of the non-linear fits are shown together 
with the error estimates calculated with the bootstrap method. 
In addition to the fit to the default data set (solid red curve), the 
dashed magenta lines show the effect of the distance uncertainty. 
Using the distance uncertainties Sd listed in Table [T] we also 
calculated the bias corrections for r(/) for distances d - Sd and 
d+Sd. Thus, the upper dashed line corresponds to distance d-Sd 
and a smaller bias correction. 

Figure [TO] shows the linear slopes k and the values of pa¬ 
rameter C for this sample of fields. The fits were performed to 
all data, without a r(/) threshold. If the data below r(/) = 0.6 
were removed, the median slope r(250yum)/r(/) = 1.6 x 10 -3 


5 Distance between a (r(/), r(250yL/m)) point and the model curve is 
measured in a coordinate system where the error region of the point 
is circular. We used the smallest distance to the curve, ignoring the 
marginal effect resulting from the curvature of the model curve. 



Fig. 10. Linear slope k (upper frame) and the parameters C 
(lower frame) for the 23 selected fields. The values obtained 
without bias corrections are shown with black symbols. The val¬ 
ues obtained with corrected r(/) and r(250yum) data are shown 
with red symbols, the shaded area corresponding to the uncer¬ 
tainty of the bias correction that is due to the uncertainty of the 
distance estimates. The dashed lines show the median values 
corresponding to the black and red symbols. The fields are ar¬ 
ranged in order of increasing distance, and the r(250/mi) values 
are based on SPIRE data alone. 


did not change appreciably (by less than O.lxlO -3 ). The me¬ 
dian value of C increased to 2.0 x 10~ 4 , which in that case is 
still lower than the scatter. The differences between the fields 
are larger than the estimated formal uncertainties (including the 
statistical errors of r(/) given by NICER and r(250yum) derived 
from the uncertainty of the surface brightness measurements). 
The error bars only reflect the statistical errors of the fits and do 
not include the uncertainty of the bias correction, for example. 
Figure [8] demonstrates that the uncertainty of the distances can 
be a significant source of error. In Fig. [TO] the shaded areas show 
the difference between the values obtained with distances d-Sd 
and d + Sd , as listed in Table [T] These were estimated directly by 
repeating the analysis using these two distance values. A smaller 
distance corresponds to smaller bias correction in r(/) and, thus, 
to a steeper slope k and typically a higher value of C. In some 
cases, the value of C obtained with the default distance d is out¬ 
side the shaded region, showing that the effect is not always this 
simple. The distance uncertainty is not yet enough to explain all 
the scatter in k and especially in C. A change in the distance es¬ 
timate results at first approximation in a nearly linear scaling of 
t(J) values (see Fig. [8] comparison of the dashed magenta lines). 
In reality, the situation may be more complex. In particular, if a 
field contains cloud structures at different distances, this might 
result in large errors in both k and C. Figure [TO] also shows that 
in spite of the small formal errors of the least-squares fits, we 
cannot constrain the opacity values in the last two fields with 
distances exceeding 2 kpc. 


In the sample of Fig. [TO] the median value of C is close to 
zero with a number of fields with negative values. The positive 
values of C in Fig. [6] are due to a small number of fields, and 
the increase of r(250yum)/r(/) values was only visible above 
r(/) ~ 5. There are only 20 fields with any data points above 
r(J) = 5. Only six fields have ten or more data points above this 
limit: G6.03+36.73, G70.10-l.69, G82.65-2.00 G92.04+3.93 
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G107.20+5.52, and G202.02+2.85. Of these, only G6.03+36.73 
is included in the sample of Fig. [K)] All the others were excluded 
because the bias correction changed the slope k by more than 
30%. Thus, a clear steepening of the relation r(250/im vs. r{J) 
is seen exclusively in fields with the highest column densities, 
which for the same reason also have the largest uncertainty re¬ 
garding the bias corrections. 


3.6. Maps of r(250yum)/r(/) ratio 

We also examined the ratios r(250yum)/r(/) in the form of maps. 
This is useful if k changes in small regions that have little effect 
when all data of a field are fitted. Unlike in Fig. [8j where the 
offset between r(250yum) and r{J) is a free parameter, the ap¬ 
pearance of the ratio maps depends on the consistency of the 
r(250yum) and r(J) zero points. Because we do not have an ab¬ 
solute zero point for r(/), we used the reference areas listed 
in Table [T] and subtracted from r(250yum) and r{J) the average 
value found in the reference area. This limits the region where a 
reliable ratio can be calculated, excluding regions of low column 
density. The details of the calculations are given in Appendix |E| 
where we also show the figures of selected fields. As an ex¬ 
ample, Fig. [TT] shows the field G4.18+35.79 (LDN 134), where 
the ratio r(250yum)/r(/) is strongly correlated with column den¬ 
sity. In Fig. |TT] the increase of k remains clear even in maps 
of (r(250yum) ± Sr{250pm) / (r{J) ± Sr{J)). The error estimates 
St( 250/rm) include the statistical errors due to Herschel pho- 
tomet ry and unc ertaint y of the surface brightness zero point (see 
Sects. [Z23 and 2.2.2). The parameter St(J) corresponds to the 
uncertainty of the r{J) zero point (see Appendix [E]). For r(J) the 
formal error estimates calculated with NICER are below 10%. 

For high-opacity sources like G4.18+35.79 the bias correc¬ 
tion of t(J) is very important. If no background stars are visible 
through some part of the core, the values of k naturally remain 
more uncertain (see, however, Sect. |3.7| ). Conversely, the zero- 
point uncertainty only becomes important in diffuse regions but 
might even reverse the correlation with column density. The ratio 
maps are also affected by the assumption of a constant value of 
P and of potential errors in the bias corrections. However, we ar¬ 
gue in Sect. |3.7| that these are mainly multiplicative errors (that 
do not affect the morphology of the maps) and/or tend to de¬ 
crease the variations seen in the ratio maps. Therefore, we are 
confident that the increase of submillimetre opacity that is seen 
in some of the maps is real. 

Based on the maps, the submillimetre opacity is correlated 
with column density in the fields G4.18+35.79, G6.03+36.73, 
Glll.41-2.95, G161.55-9.30, G151.45+3.95, and G300.86-9.00 
(see Appendix[E]). In G4.18+35.79 and G6.03+36.73 the values 
rise to close to 4 x 10 -3 . In some fields the background sub¬ 
traction reduces the available map area to such an extent that no 
conclusions can be drawn. 
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Fig. 11. Field G4.18+35.79 (LDN 134). The upper frames show 
r(250/im) (frame a) and the ratio r(250yum)/r(/) (frame b). The 
lower frames show the lower (frame c) and upper (frame d) limits 
of r(250yum)/r(/) calculated as (r(250yum)+^r(250yum))/(r(/)- 
Sr{J)) and {T{25tdpm)-5r{25f)pm)) / {t{J)+5t{J)). The areas not 
covered by Herschel observations and regions with a SN below 
0.5 have been masked. In frame a, the solid black contour and 
the dashed white contour correspond to r(250/rm = 6r( 250/rm) 
and r(/) = Sr(J). The maps have a resolution of 180 " and t(J) 
is derived using 2MASS data. 


in the extinction maps. The analysis was repeated for the four 
fields with VISTA data with resolutions of 180" and 120" . The 
results are summarised in Fig. [12] The resolution has no strong 
systematic effect on the parameters. The largest differences ap¬ 
pear when parameter C is estimated excluding low column den¬ 
sity points. However, even in that case the difference between 
the results with a resolution of 120" and 180" is smaller than 
the effect of excluding low r(J) data from the fits. When VISTA 
data are available, the results are close to those obtained with 
2MASS. Because we used the same Herschel data and bias cor¬ 
rections derived in the same way, the results are not independent. 
However, because the uncertainty of r(J) is expected to be one 
of the most significant sources of error, this gives us some con¬ 
fidence that the observed differences between the fields are real. 

The r(250yum)/r(/) ratio in the field G4.18+35.79 was 
shown in Fig. |TT] Figure 13 shows the corresponding figure 
obtained with VISTA NIR data and a spatial resolution of 
120". The highest value of r{250pm)/r(J) has increased from 
3.6 x 10 -3 in Fig. IT to 6.7 x 10 -3 . This is mainly attributed to 
the increased spatial resolution, although also at the 180" reso¬ 
lution the peak value is ~25% higher than in Fig. 11 Even in 
VISTA data there are only ~10 stars within the 2' x 2' area cen¬ 
tred on the r(250/mi maximum, and therefore the peak value of 
r{250pm)/r(J) is subject to some uncertainty. 


So far, all NIR extinction maps were calculated at 180" res¬ 
olution. Depending on the number of background stars, extinc¬ 
tion map could be derived at a higher resolution and possibly 
with smaller bias. This especially applies to the four fields for 
which VISTA observations are available. We recalculated the 
extinction maps at 120" resolution, repeating Monte Carlo sim¬ 
ulations to estimate the bias of r(/). The smaller beam increases 
the noise per resolution element, but does not yet cause holes 


3.7. Potential systematic errors 


Bec ause of the significance of the bias corrections (Sects. |3.2| 
and 3.3 ), we tried to characterise the effects that systematic er¬ 
rors in these corrections could have on the r(250yum)/r(/) ra¬ 
tios. Furthermore, the assumption of a constant dust emissivity 
spectral index may be incorrect. Below we examine the possible 
systematic effects caused by these factors. 

The bias correction made to the r(250/im) values is in it¬ 
self small (see Fig. [4]) and, consequently, the errors made in that 


10 


























M. Juvela et al.: Galactic cold cores V. Dust opacity 


-4°24'00.0" 

30’00.0" 

se’oo.o” 

42'00.0" 

48'00.0" 

54'00.0" 

-4°24'00.0" 

30'00.0" 

36'00.0" 

42'00.0" 

48'00.0" 

54'00.0" 


a G4.18+35.79 


■ 5 


14 
12 
10 

6 
4 
2 
0 

-2 








4.4 
4.0 

3.6 

3.2 
2.8 

2.4 
2.0 

1.6 

1.2 
0.8 



15h54m00s 

RA (J2000) 


15h54m00s 

RA 02000) 


Fig. 13. Ratio r(250yum)/r(/) in field G4.18+35.79 at a resolu¬ 
tion of 120", based on VISTA NIR data. Frame a shows a map 
of r(250/im), frame b a map of the ratio r(250yum)/r(/). Frames 
c and d are estimated lower and upper limits of r(250yum)/r(/) 

(cfFig.[n}. 


correction are expected to be small. The correction was derived 
from radiative transfer models with dust properties correspond¬ 
ing to R y =5.5 ( |Draine||2003b| . If the submillimetre dust emis- 
sivity is in fact higher, the dust column density will be overes¬ 
timated and models will also overestimate the cloud opacity at 
visual and NIR wavelengths, thus exhibiting stronger tempera¬ 
ture variations than the real clouds. Because the r(250yum) bias 
is related to the line-of-sight temperature variations, we could in 
this case systematically overestimate the bias in r(250yum). To 
check the magnitude of the effect, we repeated the modelling 
using a dust model with higher long-wavelength emissivity. We 
used a dust model from |Ossenkopf & Henning] ( |1994| ) with co¬ 
agulated grains with thin ice mantles accreted in 1(P years at a 
density of 10 6 cm -3 . Compared to NIR (the wavelengths con¬ 
tributing to most of the heating deep inside a cloud), this dust 
model has an emissivity higher by -50% at SPIRE wavelengths. 
The results are shown in Fig. [14] the open circles correspond¬ 
ing to this alternative modelling. Because of the smaller esti¬ 
mated bias, these should be below the k values of our previ¬ 
ous analysis (red symbols). The median value of k is 1.58xl0 -3 
and thus practically unchanged. The strongest change is seen for 
G6.03+36.73 (LDN 183), where the estimate has been reduced 
by more than 20%. This is a source with very high column den¬ 
sity and thus, in its central parts, a very uncertain estimate of 
r(250yum). However, the uncertainty of r(250yum) bias must also 
be considered in connection with the uncertainty of the r(/) cor¬ 
rection, which could compensate for some of the change (see 
below). 

The r(/) bias corrections are more significant than the 
r(250yum) bias corrections. The estimation of the r(J) bias is, 
in principle, more reliable because it only depends on the as¬ 
sumption ofthe t(J) structure of the clouds. In our calculations 
(see Sect. |3.2| ), r(/) was derived from Herschel observations, di¬ 
viding r(250/im) by the constant factor of k =1.5xl0“ 3 to obtain 
a template map of r(/). There are two possible sources of error. 
First, if the targets contain much structure below the 18" resolu¬ 
tion of the r(250yum) maps, we will underestimate the bias and 
our k estimates will be too high. We cannot directly estimate this 
error, but it is expected to be a small fraction of the total bias 
estimate. This is because 18" «: 3' and, at least for the clos¬ 
est fields, Herschel already resolves most of the cloud structure. 


The second potential source of error is again connected with dust 
emissivity. If the local ratio between r(250yum) and t(J) is higher 
than 1.5xlO -3 (strongly increased submillimetre opacity), we 
have overestimated the cloud opacity at NIR wavelengths in the 
modelling, the bias correction of r(/) is too large, and we un¬ 
derestimate the true value of k. Thus, a change in the value of k 
that is used to estimate the r(/) bias will change the recovered 
value of k in the same direction. To examine this potential prob¬ 
lem more quantitatively, we repeated the bias estimation using k 
values of l.OxlO -3 and 3.0xlCT 3 . The resulting values of k and 
C parameters are shown in Fig. [l4]as triangles. The initial as¬ 
sumption of k =1.0xl0 -3 leads to a recovered median value of 
k =1.26xl0 -3 . The initial assumption of k =3.0xl0 -3 leads to a 
recovered median value of k =1.92xl0 -3 . In both cases the input 
and output values are inconsistent, unlike in our previous analy¬ 
sis, where an assumption of 1.5xlCT 3 led to a recovered value of 
1.6xl0 -3 . Furthermore, an error in the assumed value of k leads 
to a systematic error in the recovered value that is about half of 
the original error, even lower if the true value of k was initially 
overestimated. 

The previous test shows that the estimates of k will be bi¬ 
ased towards the selected value of 1.5xl0 -3 . Our previously re¬ 
covered median value of 1.6xl0~ 3 is thus not significantly af¬ 
fected (bias lower than 0.05), but the effect can be stronger for 
individual fields. For example, in G4.18+35.79 the estimate was 
1.8xl0~ 3 , but the true value is probably higher by -10%. The 
calculations could be iterated, field by field, to carry out the bias 
correction self-consistently with the final k estimate. However, 
the errors are typically below 10%, and rough estimates of their 


magnitude can be seen in Fig. 14 


There is a specific consequence of the way the r(250/rm) 
and t(J) bias corrections are implemented. If a cloud included 
regions of such a high opacity that no 2MASS stars were vis¬ 
ible through the cloud, the ratio of r(250/im and r(J) would 
normally be overestimated. The resulting apparent increase of 
submillimetre opacity could thus be an artefact resulting from 
errors in extinction values. However, in our analysis we also cor¬ 
rect t(J) in this case based on the assumed opacity derived using 
the r(250yum) input map and the extinction calculated with sim¬ 
ulated 2MASS stars. If the gap in the distribution of background 
stars is increased, the ratio r(250/im)/r(/) does not continue to 
increase, but instead tends towards the assumed ratio, 1.5xl0 -3 . 
Higher values should thus not be the result of gaps in extinction 
data. 

To investigate the potential effects of a spatially varying 
spectral index, we repeated the analysis using an ad hoc J3(T ) law 
to introduce p variations in all of our maps. We took the temper¬ 
atures calculated with p = 2.0 and fixed new y 3 values pixel by 
pixel using a functional dependence y 3 = 2.0 x (7715.0) -0 - 24 . We 
then repeated the full analysis, starting with the zero point and 
colour corrections and continuing with the calculation of colour 
temperatures and submillimetre opacity. We did not solve the 
(T, p) values (which are very susceptible to noise effects), but 
simply assumed that p could vary in a systematic way so that 
the values are higher when the dust temperature is lower. The 
parameters of the fi(T) formula were selected so that ft changes 
from - 1.8 in warm regions with T ~ 23 K to - 2.2 in the cold- 
est spots T ~ 10 K (cf|Dupac et al.||2003[ |Desert et al.||2008 


Planck Collaboration XXIII 201 1 HParadis et al.|2010HVeneziani| 
et al.|2010[ Juvela et al.|201 1| ). The average p is still close to the 
original p = 2.0, and we mainly examined the effects of corre¬ 
lated changes of p rather than the effects of absolute p values that 
can be estimated more directly. The crosses in Fig. 14 show the 
slopes k = Ar(250/im/Ar(/) and the parameters C obtained with 
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Fig. 14. As Fig. [KM but comparing our default r(250yum)/r(/) 


estimates (red solid circles) and results from alternative anal¬ 
yses: r(250yum estimates derived assuming a spatially varying 
spectral index (cro sses), r(250yum bias estimated with Ossenkopf 
& Henning (1994) dust model (open circles), r(/) bias estimated 
using k- 1.0 x 10 -3 (triangles pointing upwards), and r(/) bias 
estimated using k- 3.0 x 10 -3 (triangles pointing downwards). 


the spatially varying fi. The values are typically within -10% of 
the values obtained with J3 = 2.0. The strongest changes of k are 
seen in the two closest fields, G4.18+35.79 and G6.03+3673, 
which have well-resolved clumps with very low temperatures. 
The increase of the slope values is -20%. Note that this is mostly 
consistent with the general dependence between r(250yum) and 
and not necessarily an effect of the spatial variation of fi. In this 
respect, it is very interesting to note that the effect on the param¬ 
eter C is weak. In other words, if the variations of are as as¬ 
sumed above, this will be reflected in the slope r(250/im)/T(/), 
but without a noticeable non-linearity in the r(250yum) vs. r(/) 
relation. 

All the above suggests an uncertainty of 10-20% in k and 
-0.1 units in C. In particular, if the increase of dust opacity is 
associated with values /3 > 2.0, our highest k estimates could 
still systematically underestimate the true values of k because of 
the lower assumed value of and because the r(/) correction 
biases the k values towards 1.5 xlO -3 . 


4. Discussion 

We have examined the submillimetre opacity by correlating the 
250 yum optical depth r(250/im) with the near-infrared optical 
depth, r(/), assuming the latter to be an independent tracer of 
the total dust column density. Because the comparison was made 
at a resolution of 2' , it is not sensitive to dense cores (Ay » 
10 mag), at which both the r(250//m) and r(/) estimates would 
become very uncertain. Nevertheless, corrections for systematic 
bias in r(250/im) and especially in r(/) are important. 

The sample consists of the heterogeneous set of 116 Galactic 
fields that were mapped with Herschel as part of the Galactic 
Cold Cores project. The main objectives were to estimate the 
typical ratio of r(250/im)/T(J) and to search for variations of this 
quantity, between the fields and as function of column density. 
Such variations were then related to differences in the properties 
of interstellar dust grains. For the present sample, high column 
densities also imply low dust temperatures and thus conditions 


where submillimetre dust opacity is expected to be enhanced by 
grain aggregation. The limited resolution means that we did not 
probe the full range of opacity variations, if these are partly lim¬ 
ited inside compact cores. 


4.1. Main results and their reliability 

By restricting the analysis to -20 fields for which the re¬ 
sults appeared most reliable, we derived a m edian value of 


r(250yum)/r(J) = 1.6 x 10“ 3 (see Sect. 3.5 and Fig. 10). In 


Fig. [10] 50% of the most nearby fields had positive values of C, 
the multiplier of the second-order term, indicating some degree 
of positive correlation between column density and submillime¬ 
tre emissivity measured by r(250yum). In the maps of the ratio 
k , the same tendency was very clear in only six cases. The low 
percentage is partly caused by the noise in r(/), whose mag¬ 
nitude is strongly correlated with the cloud distances. For two 
of the best examples, G4.18+35.79 and G6.03+36.73, the ratio 
r(250/im)/T(J) increases to - 4 x 10 -3 , almost a factor of three 
higher than the median value. The peak values are uncertain be¬ 
cause of the large bias corrections, but because of the low spatial 
resolution used, the strongest effect can be even greater. 

No dependence on either Galactocentric distance or on 
Galactic height was observed (Fig. [5]). The reliability of the es¬ 
timates decreases with distance, but nevertheless, our sample 
extends over more than ~4kpc in Galactocentric distance. No 
trends are seen; if they were present at 10% level, they should 
still be visible over the scatter of individual data points. The re¬ 
sult might be affected by systematic errors in the bias correction. 
However, these probably depend either on the distance or on the 
morphology of the field and at first approximation are not ex¬ 
pected to be different in the inner and in the outer Galaxy. 

In Fig. [lO] the only clear trend is the decrease of k as a 
function of distance (also visible as larger scatter around the 
Gal acto centric distance of 8.5 kpc). As mentioned in Sects. 3.2 
and |3.3| there are two possible explanations. First, the bias cor¬ 
rection of t(J) might be overestimated so that as the distance in¬ 
creases, the error increases and k becomes underestimated. This 
is probably not the main reason because the bias correction de¬ 
pends as much on column density values and column density 
gradients as on distance. The trend probably is a combination 
of selection and resolution effects. At 1 kpc the 180" resolution 
corresponds to almost 1 pc in linear scale. Therefore, we mea¬ 
sure the mean cloud properties for the most distant fields. In the 
nearby fields individual clumps are resolved and the slopes k 
reflect more the contrast between diffuse regions and compact, 
core-sized objects. Thus, the trend would be compatible with the 
hypothesis that r(250yum)/r(/) increases in the densest and cold¬ 
est regions of interstellar clouds. 

The global non-linear fits of Fig. 6] also indicated an increase 
of the ratio k = r(250yum)/r(/) as the function of column den¬ 
sity. The plots show clear deviations from a linear dependence 
beyond r(/) - 4. For r(/) - 10, k is already twice as high at 
low column densities. The results are dependent on a few fields 
with the highest column densities for which the bias corrections 
reach about ten per cent. However, if the distance dependence of 
Fig. [5] means that the bias correction of r(/) is probably over¬ 
estimated and not underestimated, the true values of k might be 
even higher. 

Figure [TO] concentrated on selected fields for which linear 
and non-linear fits were more reliable than on average. Within 
this subset, no clear distance dependence was visible in k val¬ 
ues, but the scatter is larger than the estimated uncertainties. The 
linear slopes are sensitive to the highest column densities within 
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a field. The parameter C of the non-linear fits is expected to be 
even more sensitive to these data points, which probe the varia¬ 
tions of r(250yum)/r(/) inside the fields. For the whole sample, 
the median value of C was very close to zero. Nevertheless, it 
may be significant that when the two fields at kiloparsec dis¬ 
tances are excluded, the value of C appears to decrease some¬ 
what systematically with the distance of the field. There is some 
preference for the most nearby fields to have positive values; 
here, the densest parts of the clouds are resolved. 

Two of the first fields with high k values are G6.03+36.73 
and G4.18+35.79, better known as LDN 183 and LDN 134. At 
the Herschel scale, the clouds have a simple morphology, each 
consisting of a single column density maximum that is partially 
reso lved by the 180" beam. This is illustrated by Figs. m and 
E.l The ratio k - r(250pm)/r(J) closely follows the morphol¬ 
ogy of the column density distribution, making these the best 
examples of clumps with increased submillimetre opacity. The 
highest values are ~ 4 x 10 -3 , almost three times the average 
value over all fields. 

In our results, some of the main sources of uncertainty are 
the corrections made for the expected bias in r(/) and, to lesser 
extent, in r(250yum). According to Fig. [10| the net effect of all 
bias corrections is a -20 % decrease in the value of k. This num¬ 
ber applies to nearby fields but is more dramatic if all fields 
are taken into account. Figure [5] shows that for fields at - 1 kpc 
distance the correction is almost a factor of two. The correc¬ 
tions have been remarkably successful in decreasing the scatter 
of r(250yum)/r(/) values, which strongly suggests that they are 
approximately of the correct magnitude. 

In the statistical sense, estimating the r(J) bias is straight¬ 
forward, using a model of the Galactic stellar distribution and 
the higher resolution Herschel data as a template for t he co lumn 
density structure of the field (see Appendix [B] and Sect. |3.7| ). The 
correction is large, but on average, the correction itself probably 
does not suffer from major systematic errors. Errors could arise 
either from incorrect distance estimates or from the use of an in¬ 
correct model of the NIR opacity distribution. The distances are 
uncertain, and through the r(J) bias, their effect on the main pa¬ 
rameters is shown in Fig.[K)](the grey bands) and in Fig. [8] For a 
sample of fields, this is mainly a statistical and not a systematic 
error. 

The model of r(J) distribution in each field was derived from 
Herschel data at 18" resolution. If there is still significant struc¬ 
ture below this scale, our correction of t(J) values will be too 
low. The difference of the 2-3' scale and the 18" scale is so large, 
however, that most of the effects of column density variations 
are already included. Thus, after the distance, the other main er¬ 
ror in the r(J) bias correction arises from scaling the Herschel 
estimates of r(250yum) into a template of the J-band opacity. 
We showed in Sect. T7 that this amounts to an uncertainty of 
-10% in the final listed values of k - Ar(250yum)/Ar(/). For 
the sample in Fig. 10 the value of k assumed during the bias 
correction is consistent with the recovered value (to be precise, 
the assumed value of k = 1.5 x 10 -3 , the recovered value of 
k = 1.6 x 10 -3 ). For individual fields, the values are biased to¬ 
wards the initial assumed value. Figure [l4| showed that if bias 
correction assumed a value of k = 1.0 x 10 , the recovered me¬ 
dian value was still higher than k - 1.2 x 10~ 3 . This shows that 
k cannot be significantly overestimated because of an erroneous 
bias correction in r(/). Thus, the result that the average ratio 
k = Ar(250yura)/Ar(/) is clearly higher than the normal value 
found in diffuse medium remains robust. 

Because of the bias in the r(/) correction, our estimates are 
probably conservative for fields where values k > 1.5xlO -3 were 


obtained. The dependence on the assumed value of k decreases 
as the true value of k increases. This is because higher k means a 
lower NIR opacity and lower overall bias in r(/). Nevertheless, 
for fields like G4.18+35.79, we might systematically underesti¬ 
mate k by -10% because of this error in the r(J) bias correction. 
At the highest column densities Herschel emission data them¬ 
selves will underestimate the clo ud opacity, which leads to the 
corrections discussed in Sect. |3.3| For the optical depth ratio, the 
effect of r(250yum) bias is weaker than that of r(/). Nevertheless, 
the errors made in the corrections of r(250/mi) and r(/) (for ex¬ 
ample, those associated with a change of submillimetre opacity) 
would partly cancel each other out. 

In the analysis we assumed that apart from the problems as¬ 
sociated with the sampling provided by the background stars, 
NIR reddening is an independent and reliable measure of col¬ 
umn density. Unlike in the optical range, the NIR extinction 
curve is often assumed to be constant for a wide range of col¬ 
umn densiti es ( [Cardelli et al.|1989}|Martin & Whittet|1990||Roy 


|et al.||20l3] ). Nevertheless, some cloud-to-cloud variations are 
observed, the ratio Ej- H IE H - K ranging from values lo wer than 
1.5 to higher than 2.0 ( [Racca et al.[[2002[ |Draine|2003a ). Clear 
changes take place at high optical depths, above Ay - 20 mag, 
but only in the form of the flattening of the MIR extinction curve, 
at wavelen gths above 3 jjm pndebetouw et al.|2003}|C ambresy 
|et al.|20lH|Ascenso et al.|2013| ). Recently, |Whittet et al.| ( |2013| ) 

studied cloud LDN 183 (which is also included in our sam¬ 
ple). Comparison with the 9.7 /im silicate absorption feature sug¬ 
gested that the NIR colour excess might not be a perfect tracer 
of the total dust column. The ratio between Ej- K and the 9.7 yum 
feature was observed to increase as the column density exceeded 
Ay -20 mag. The observation might partly arise because the 
9.7 //m feature is dampened by the formation of ice mantles (see 
also [Chiar et al.|2007| . However, if we assume that NIR extinc¬ 
tion indeed overestimates the total dust column in this object, the 
ratio of r(250 /im) relative to column density would increase by 
-20% (see |Whittet et al.[|2013[ Fig. 12). This is still a weaker 
effect than the increase by more than a factor of two that was 
observed in r(250yum)/r(/). We recall that at high column densi¬ 
ties, above Ay - 20 mag, the values of r(250yum) are also uncer¬ 
tain and can be underestimated by a significant fraction ([Pagani 


et al.| ( |2004| ) discussed a similar limitation at the nearby wave¬ 


length of 200yum). 

Finally, we recall that r(250yum) values were derived us¬ 
ing a fixed value of the spectral index, = 2.0. By assuming 
P = 1.8 instead, the r(250yum) and k values would decrease 
by up to 30%. Conversel y, if the value of increased towards 
dense and cold regions ([Pupae et al.||2003[ Desert et al.||2 008 


Planck Collaboration XXIII 2011 [ |Paradis et al. |20 1 0[ | Veneziani| 
|et al.|2010[ | Juvela et al.|2011| ), we underestimate the dust opac¬ 
ity changes if we use a constant value of (3 . In this sense (and 
regarding possible bias in r(/) corrections), Figs. [6]-8 give con¬ 
servative estimates of the possible increase of submillimetre dust 
opacity. Clearly, the assumed values of y 6 must also be taken into 
acco unt when comparing our results with other studies (see Sect. 
|4.2| ). On the other hand, our tests indicate that spatial variations 
of p probably do not have a strong additional effect on k (see 
Sect. [3771 


4.2. Comparison with other studies 

The submillimetre dust opacity has previously been studied in 


et al. 2009; Juvela et al. 2011; Martin et al. 2012; Roy et al. 

2013J Malinen et al. 

2013 

2014) and to HI (Boulanger et al. 
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1996; Lagache et al. 1999 

Planck Collaboration XXIV 

2011; 

Planck Collaboration XXV 

2011J[Martin et al. 2012). Bou 

anger 

et al. (1996 ) compared COBE observations of FIR dust emission 


ered relatively constant in the NIR regime (e.g . Indebet ouw et al 
20051 Lombardi et al.|2006] Roman-Zuniga et al.|2007[|Ascenso 


with the hydrogen 21 cm observations. At high latitudes, where 
HI is a good tracer of the full gas column density, the derived 
value was t/N h = 1.0 x 10" 25 (A/ 250/im) cm 2 H _1 . Similar val¬ 


et al. 


the N: 


2013[ Wang & Jiang 2014) and the power-law slope of 
R extinction curve typically only varies at a level of -5% 


( [Stead & Hoare|2009[ [Fritz et al.|2011| . In very dense cores the 


ues were obtained in the first studies using Planck data (Planck 
[Collaboration XXIV|201 j) . 

In the most recent Planck papers, somewhat lower 


formation of ic e mantles and the grain growth could have an ad¬ 
ditional impact. Roman-Zuniga et al. (2007) examined the cloud 


Barnard 59 up to Ay = 59 mag, but found no significant changes 
in the NIR extinction law (compatible with Ry = 5.5). Similarly, 


0.55 x 10- 25 cm 2 H 1 ( 

Planck Collaboration XI 

2014, 

Planck the reddening NIR law up to E(H- K) - 1.5 mag (Ay - 9 mag). 

[Collaboration Int. XVI 

; 2014). The change is associated with Therefore, it is not likely that our results are significantly biased 


the revised calibration of the highest frequency channels and, 
correspondingly, a lower value of /3. In [Planck Collaboration] 
|Int. XVII] ( |2014| ) the spectral index above 353 GHz was found 
to be p - 1.65, lowe r tha n the value of 1.8 assumed in 
|Planck C ollaboration XXIV] ( |2011| ) or the value of 2.0 used 


of our extinction maps, practically all our data points are lower 
thanAy - 10 mag and are thus not affected by extreme optical 
depths. However, Whittet et al.| ( |2013| ) observed a decrease in 


in 


Boulange r et al.| (|1996| ). The higher dust opacity value of 


Boulanger et al. (1996) is thus largely explained by the higher 


value of 13 . 

The Planck result for the diffuse medium can be compared to 
our result most directly by converting their Ah to r(/). We can 
use the conversion factor A(H 2 )/A V = 9.4x 10 20 cm -2 mag -1 de- 
rive d by|Bohlin et ak|(|1978|) at l ow extinction, E( B- U)j0.5 (see 
also Nozawa & Fukugita 2 013|). Some studi es ( Rachford et al 


2009[ |Planck Collaboration xT]|2014[ |Liszt||2014| ) have found 
Nu/E( B - V) values that are 10-30% higher than in |Bohlin et al.| 
( |1978| ), these results apply partly to even more diffu se lines of 
sight ( E(B - V) ^ 0.1). On the other hand, |Gudennavar et al.| 
( |2012| ) examined a sample of lines of sight with E(B - V ) ex¬ 
tending to values higher than one. The result, N(H)/E(B - V) = 
(6.1 ± 0.2) x 10 21 Hcm~ 2 m ag -1 , wa s close to that of Bohlin 


|et~aL| (1978]). With the IBohlin et al. 


3.1 extinction curve (Cardelli et al. 
T(250^m)/A H ~ 0.55 x 10“ 2:) cm 2 H- ] 
(2014) corresponds to r(250/im)/ r{J) 


1978) relation and Ry = 
1989), the Planck result 


Planck Collaboration XI 
0.41xl0 _ L Our median 
1.6x 10 -3 is thus 3.9 times higher, and 

3 


value of r(250/im)/r(/) 
our highest local values, 4.5xlO -3 in G4.18+35.79, are higher 
by one order of magnitude. In [Planck Collaboration XI[ ( |2014| ) 


the estimated value of N(R/E(B - V )) was -20% higher than 
in |Bohlin et al.| ( |1978| ). With this value, our median value of 
T(250yum)/r(/) would be - 3.3 times higher than the Planck es¬ 
timate. 

We can alternatively convert our result into r(250yum)/An 
, but in dense regions the shape of the extinction curve (de¬ 
pendence on Ry) and the ratio between visual extinction and 
total column density are more uncertain. Using the Ry = 3.1 


LI83 in the ratio of E(J - K) and the 9.7 /a m silicate absorption 
feature. Above E(J - K) - 1.0 (Aj - 1.7 mag, Ay - 5.0mag 
for Ry = 5.0), the ratio deviated from diffuse medium value by 
- 20%. Stronger deviations were only seen beyond E(J-K) - 3 
(Aj - 5 mag, Ay - 15 mag for Ry = 5.0). This is above the 
range probed by our measurements. It is not clear that changes 
in this ratio would only be caused by the NIR extinction curve. 
However, if, in some sources, the extinction curve does flatten 
above A(/) - 2, this might contribute to the observed increase 
of r(250yum)/r(/) (possibly at the 20% level). 

The effect of the bias corrections on r(250yum) /t(/) is typ¬ 
ically -20% or weaker, and their uncertainty is smaller. The 
largest uncertainties in r(250jum)/AH are caused by /3 and pos¬ 
sibly by Ry. With /3 = 1.8 and R v =5.5, the r(250/im)/r(/) 
values would be 40% lower than with /3 = 2.0 and Ry=3.1. 
At this lower limit, our median value of r(250yum)/r(/) would 
not be 3.9 times, but -2.3 times the value found in the diffuse 
high-latitude sky. The 40% uncertainty may be a realistic 1 cr 
lower limit for the linear fits that include all pixels in the maps. 
However, it is a conservative estimate for the clumps where the 
average value of /3 is expected to be clearly higher than 1.8. In 
fact, preliminary results indicate that the average value of /3 in 
our fields (including the more diffuse regions) is close to/3 = 1.9, 
and that with f3 = 2.0 we underestimate t he r(250yum) values of 
many clumps ( [Juvela|2014 in preparation) ). 


Increased far-infrared and submillimetre opacity has been 
reported by many authors (Kramer et al.||2003[ Lehtinen et al 


2004; del Burgo & Laureijs 2005; Ridderstad & Juvela 2010; 


Bernard et al.||2010[ |Suutarinen et al.||2013[ etc.). One famous 

example is the Taurus filament LI506, for which the models 
of |Stepnik et al.[ ( ]2003| ) suggested an increase by more than a 


corresponds to r(250/im)/AH = 2.16 x 10 25 cm 2 H 1 (again, 

properties and the structure of this filament, |Ysard et al. 

|2013) 

of course, 3.9 times the Planck value). However, Ry is ex- 

estimated the increase of the 250 jam opacity to be -2. 

Martin 


pected to be higher than 3.1 in dense clouds. Using Ry = 5.5 
instead of Ry = 3.1 in converting r(J) into visual extinc¬ 
tion, the values of r(250jum)/AH would decrease by - 15% 
(change in Ay/E(J - K )). However, the scaling we used above 
(corresponding to the Ry=3.1 extinction_curve and the ratio 
of A(H 2 )/A V taken from Bohlin et al.| ( [1978 )) corresponds to 


N(H)/E(J - K) = 11.0 x 10 Z1 cm z mag L which is very close 
to the value of 11.5 x 10 21 cm -2 mag -1 that Martin et al. ( 2012 ) 


derived in Vela molecular cloud using 2MASS NIR data up to 
E(J - K) - 0.55 mag (E(B - V) - 1.1 mag). 

When Ry is modified, the NIR extinction curve remains 
practically unchanged and the differences take place between 
the optical and NIR wavelengths. Compared to shorter (optical) 
and longer (MIR) wavelengths, the extinction curve is consid- 


BLAST and IRAS data with the reddening of 2MASS stars. 
The properties of the examined areas corresponded to the av¬ 
erage properties of our fields, with column densities extending 
to 10 22 cm -2 and with typical dust temperatures of - 15 K. They 
found a very similar range of dust opacities, T(250/dm)/Nu = 
(2 - 4) x 10 _25 cm 2 H -1 (assuming /3 = 1.8). In the Orion A 
cloud, a comparison of Herschel and 2MASS data led to the de¬ 
tection of a dependence of A 028 on the 250yum opacity (Roy 
et al.|2013| . The range of column densities in Orion A was sim¬ 


ilar to our fields, and the derived dust opacities were mainly in 
the range of r(250yum)/AH = (1 - 3) x 10 -25 cm 2 H -1 . These es¬ 
timates were derived with (3 = 1.8 and would become - 20% 
higher (depending on the details of the fitting) if a value of 
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P = 2.0 were used. More recently, Lombardi et al. ( |2014| ) de¬ 
rived ratios A(K)/r(850yura) of 2640 mag and 3460 mag for the 
Orion A and B molecular clouds, respectively. The 850/im opti¬ 
cal depth was derived from Herschel observations rescaled using 
a comparison with Planck measurements, and the NIR extinc¬ 
tion was calculated with the method NICEST ( |Lombardi[2009| . 
T he modified blackbody fits used ft values that were estimated 
in [Planck Collaboration XI1 ( [2014| at a resolution of 30'. With 
the reported average value of ft = 1.8 and adopting a ratio 0.40 
between A(K) and A(J), the results correspond to r(250yum)/r(J) 
values of 1.56 x 10 -3 and 1.19 x 10 -3 for Orion A and B, re¬ 
spectively. These fits were made for data r(850yum) < 2 x 10~ 4 
(r(250yum) ^ 1.8 x 10 -3 ). The optical depth range is similar to 
many of our fields, and the result for Orion A is close to our 
median value for the fits concerning entire fields. By adopting 
P - 1.8, our median value would fall between the Orion A and 
Orion B estimates of Lombardi et al. ( 2014]). 


4.3. Implications for dust evolution 

As shown by simulations, an observed increase of the dust far- 
IR/submm opacity towards dense regions cannot be due to ra¬ 
diative transfer eff ects, but must originate in intrinsic variations 
of dust properties ( [Malinen et al.||20lH | Juvela & Y sard|[20l2] 
|Ysard et al.|2012| ). Theoretical studies have shown that such an 
increase, coupled with a decrease in dust temperature , can be ex- 
plained by the formation of large aggregate particles ([O ssenkopf 
& Henning||1994[ |Stognienko et al.||1995[ |Ormel et al.||2011| 


Koh ler et ai.|20li[|2012| ). |Kohler et al.| ( |2bi2| ) showed that the 

coagulation of just four big grains of the diffuse ISM type, coated 
by smaller carbon grains, already leads to an increase in the 
opacity at 250//m of a factor 2.6. These authors also showed 
that these aggregates can form within a typical cloud lifetime of 
10 million years ( [ Walmsley|| 19911 ). Consequently, we interpret 
the observed increase in r(250yum) in our cold core sample as 
grain growth in dense molecular regions. 

In our sample, the clouds LDN 183 and LDN 134 
(G6.03+36.73 and G4.18+35.79) were the most convincing ex¬ 
amples of increased submillimetre dust opacity. LDN 183 has 
been studied thoroughly in both continuum and line emission 
(e.g. |Juvela et al.||2002[ |Pagani et al.||2005| |2007| ). A slight 
increase of 200pm op acity was already reported based on 
ISOPHOT observations ( [Juvela et al.|2002| ), but ISOPHOT data 
A < 200pm and, to some extent even Herschel observations, 
are not sufficient to fully probe the inner parts of the cloud 
where the high visual extinction is approaching_L00 mag _and 
the dust temperature drops well below 10 K (Pag ani et al.|2 004, 
2014, submitted). LDN 183 was the first object where enhanced 


mid-infrared (MIR) light scattering, the s o-called coreshine phe- 
nomenon, was detected in Spitzer data (Steinacker et al. 2 010| 
Pag ani et al. 2010| ). The effect was also seen in | Juvela et al.| 
( |2012| ), where WISE MIR observations were analysed and both 
LDN 183 and LDN 134 were found to be sources of strong MIR 
emission. If the signal is interpreted as scattering of the interstel¬ 
lar radiation field, it seems to imply the presence of very large, 
micrometre-sized dust particles ( Steina cker et al.||2010| |2014| 
|Lefevre et al.||2014] ). Thus, our detection of increased submil¬ 
limetre opacity in these clouds agrees with the evidence provided 
by MIR wavelengths. 


5. Conclusions 

We have examined dust optical depths by comparing measure¬ 
ments of submillimetre dust emission and the reddening of the 


light of background stars in the near-infrared. The goal was to 
measure the value of dust submillimetre opacity and to search 
for variations that might be correlated with the physical state 
and the environment of the cloud. The study led to the following 
conclusions: 

- For a subsample of 23 fields with well-defined correlation 
between the two variables, we obtained a median ratio of 
r(250yum)/r(/) = (1.6 ± 0.2) x 10" 3 . This is more than three 
times the value that was derived from Planck data for the 
diffuse medium at high Galactic latitudes. Assuming p = 1.8 
instead of p - 2.0, the value decreases by 30%, but is still 
more than twice the diffuse value. 

- The conversion to r(250pm)/Nu involves more assumptions. 
Using the Ry = 3.1 extinction curve and the A/ITU)/Ay ratio 
of|Bohlin et al.|(|l978|, our median estimate corresponds to 
r(250pm)/N H = 2.16 x 10" 25 cm 2 H -1 . 

- The fit to all data above r(J) = 1 gives a relation r(250pm) ~ 
1.25xl0 -3 r(/)+0.11xl0 -3 r(/) 2 . The positive second-order 
coefficient C = 0.11 x 10~ 3 is determined by a small number 
of fields that, because of the high column density, are subject 
to large uncertainty in the bias corrections. 

- For the same sample, the scatter in the coefficients of 
the second-order terms C is ~ 2 x 10 -4 and the median 
value is consistent with zero. Spatial variations of the ratio 
r(250pm)/r(J) are only seen in a few fields. 

- From the maps of r(250yum)/r(/), we identified six fields 
where the ratio appears to increase further at the location of 
the main column density peaks. The highest values in the 
fields G4.18+35.79 and G6.03+36.73 are ~ 4 x 1(T 3 at the 
resolution of 180". Thus, although the densest clumps are 
associated with the largest uncertainties, we consider this in¬ 
crease of submillimetre opacity to be real. 
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G6.03+36.73, d = 110pc / 


G4.18+35.79, d = 110pc 


G21.26+12.11, d = 120pc 
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Fig. 8. Fits of r(250yum) vs. r(/) in selected fields ordered by increasing distance. The red and blue points (dust temperature above 
and below 14 K) with error bars are the bias-corrected data points, where r(250yum) is based on SPIRE data. The slopes of linear 
fits are listed in the upper left corner for (1) the default data set based on SPIRE data alone (“Def”), (2) 160-500yum data (“T=4”), 
(3) SPIRE data but without bias corrections (“-deb”). The linear fit of the default case is shown with a black line. The non-linear 
fits are shown with solid blue lines (MCMC) and solid red lines (bootstrapping) with associated shaded 68 % confidence regions. 
The dashed magenta lines correspond to different bias correction of r(/) using distances d - Sd and d + Sd. The parameters from 
bootstrapping are given in the lower right corner. The non-linear fit to data without bias corrections is plotted with a solid green 
curve (without error region) with the parameter C given at the bottom of the figure. The zero points of the r(J) axes are not absolute. 
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Fig. 9. Continued... 
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Fig. 12. Comparison of parameters k (upper 
frames) and C (lower frames) obtained with 
data at resolutions of 180" and 120" . Frames 
a and c show fits to all data, frames b and d 
fits to values r(/) >0.6 alone. The open circles 
and crosses show the estimates obtained with 
2MASS data at resolutions of 180" and 120" . 
The filled triangles show the results for VISTA 
observations, the larger triangles corresponding 
to a resolution of 180" , the smaller to a reso¬ 
lution of 120" . The fields are the same as in 
Fig. [TO] with the addition of G358.96+36.75, 
which has fewer data points above tj > 0.6 and 
for which parameter C could not be fitted with 
extinction maps with a resolution of 120" . 
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Notes. The aliases refer to t he follow ing catalogues: Arch = |Dcsert et al. |([2008), B = |Barnard|([1919];|Barnard et al.|{1927], BDN = |Bernes |(|1977], CB =|Clemens & Barvainis|(|1988], FeSr = 
|Feitzinger & Stuewe|(|1984] HMST Hartley et al.|(T986] , LDN = |Ly nds | < [1962J , MBM = |Magnani et al.| ( |1985], MLB^Myers et al.| ( ll983^SLDN =|Sandqvist & LindroosH1976] . Most entries 
are looked up from |Dutra & Bica (2002]. Names starting with PCC and the entries S1-S10 refer to names used in |Juvela et al . ( 2010^ and Planck Collaboration XXII ( 201 lj for some of the Herschel 
fields. 
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Appendix A: Alternative data analysis 


To investigate the robustness of the results to details of the data 
reduction, we calculated a set of alternative r(250/im) maps. In 
addition to the analysis of three or four Herschel bands (see 
Sect. |2.2| ), we considered the use of local background subtrac¬ 
tion and the possibility of making a correction for residual map¬ 
ping artefacts with the help of other all-sky surveys. The result¬ 
ing r (25Q /jm)/r(/) ratios were compared to the values found in 
Sect. |3.1| The comparison was carried out without bias correc¬ 


tions, comparing the results with the 160-500 pm fits of Sect. 2.4 


At. Local background subtraction 

Our default analysis is based on Herschel data for which the ab¬ 
solute zero points were derived from a comparison with Planck 
and IRAS maps. As an alternative, we used Herschel surface 
brightness maps from which the diffuse background was sub¬ 
tracted using the reference regions listed in Table[l] Thus, the av¬ 
erage surface brightness of the reference region was subtracted 
from each Herschel surface brightness map separately, before 
calculating the colour temperatures and the dust optical depths. 

The local background subtraction might be a more reliable 
way to ensure a consistent zero level for the compared quanti¬ 
ties. However, it also means that colour temperature and column 
density can only be estimated in the part of the map in which 
the surface brightness values are significantly higher than those 
of the reference area. We masked the area in which the signal 
is lower than twice the estimated statistical uncertainty of the 
surface brightness in each band. The final mask is a combina¬ 
tion of these masks and the original mask that eliminated the 
map boundaries for which the convolution to the resolution of 
the t(J) data is only poorly defined. 


A.2. Checks for mapping artefacts 

Although Herschel data are usually of very good quality, there 
can still be some small artefacts that affect some parts of the 
maps. Errors might result from data reduction or from instru- 
mental effects such as strip ing or general gain changes ( |Xu et al. 
2014; Paladini et al.|2013| ). If processing includes high-pass fil¬ 
tering, the large-scale surface brightness gradients may be af¬ 
fected and the contrast between faint and bright regions may be 
decreased. Our maps often contain significant emission up to the 
map boundary. Without a flat border region with very low emis¬ 
sion, it is difficult to estimate whether the baseline assumed for 
the scans is correct. Such effects could be more important and 
more difficult to detect for small maps. Thus, this could mostly 
affect PACS maps, for which the signal-to-noise ratio is also typ¬ 
ically lower than in the SPIRE data (J uvel a et al. 2010| ). 

These effects were investigated with the help of independent 
FIR and submillimetre data. At 100 pm, 350 pm, and 500 pm we 
can compare Herschel data almost directly with the correspond¬ 
ing IRIS and Planck bands. The Planck 545 GHz data were cor¬ 
rected to 500 pm using a modified black body with the Herschel 
colour temperature map and /3 equal to 2.0. Keeping the other 
parameter constant, an error of 2 K in temperature or an error of 
0.2 in p would both correspond to only a ~2% error in the extrap¬ 
olated value. At 160 pm and 250 pm we used values interpolated 
from IRIS 100 pm, AKARI140 pm ( jMurakami et al.|2007| ), and 
Planck 350 yum channels. Here A ft ~ 0.2 translates into a change 
of less than 1% at 160yum. For a direct extrapolation from the 
350 pm to 160//m, an error of A/3=0.1 would result in an error 
in the interpolated value that is still lower than 8%. Typically 


Table A.l. Comparison of the mean values of r(250yum)/r(/) 
obtained with different versions of Herschel data, without bias 
corrections. The columns are (1) data version, (2) number of 
fields with < 10% error in k, (3) mean and standard deviation 
for that sample of fields, (4) difference and standard deviation 
when compared to the values in the default case (3-bands), for a 
common set of 81 fields). 


Data set 

Fields 

r(250yum)/r(/) 

(io- 4 ) 

A(T(250jum)/r(7)) 

(io- 4 ) 

3-band 

105 

22.4+9.6 

- 

4-band 

83 

22.2+12.5 

-0.04 + 2.80 

Bg-subtracted 

102 

20.1+9.5 

-0.96 + 1.88 

3-band, corr. 1 

105 

20.7+8.6 

-0.91 + 1.51 

4-band, corr. 1 

82 

21.3+11.8 

-0.54 + 2.43 


1 Correction at large scales using ancillary data. 


interpolation errors should thus be below the statistical errors. 
The same considerations apply to the zero-point corrections of 
Sect 2.2.2| with the difference that they are not affected by any 
multiplicative errors. 

The reference data were compared with the original Herschel 
maps at 6' resolution to derive an additive correction that leaves 
the median value of the maps unchanged and only affects scales 
larger than ~ 6'. We also checked similar multiplicative correc¬ 
tions, assuming that the zero points of the different surveys are 
compatible, and even calculating corrections where linear fits be¬ 
tween Herschel and reference data were first used to estimate the 
differences in zero-point and gain calibration. The last alterna¬ 
tive would avoid the assumptions of consistent calibration and 
zero points between the data sets. In most cases, there are no 
significant differences between the three choices. A typical map 
has no clear artefacts, and the proposed correction will prob¬ 
ably decrease the data quality. However, when the local arte¬ 
facts are clear (e.g., excessive surface brightness some corner of 
a Herschel map), the corrected map should give a better descrip¬ 
tion of the true surface brightness. Thus, we do not believe that 
the corrected maps represent a clear improvement, but the dif¬ 
ference between the corrected and uncorrected data should give 
some idea of the potential effect that artefacts of that magnitude 
could have. 


A.3. Comparison of the data sets 

We compared the r(250pm)/r(J) values (without bias correc¬ 
tions) among six data sets: (1) three bands at 250-500 pm (our 
default data set), (2) four bands at 160-500 pm, (3) three bands 
with local backg roun d subtraction, (4) three bands with the cor¬ 
rections^ Sect. |A.2[ and (5) four bands with the corrections of 
Sect. |A.2| Table | A.l] shows the results, comparing the disper¬ 
sion of the obtained r(250yum)/r(/) values between fields and 
the change in the values compared to the default case where three 
bands and the absolute zero points were used. 

The second column of the table lists the number of fields 
where the formal error of the slope r(250yum)/r(/) is lower than 
10%. This number is smaller when PACS data are included, 82- 
83 fields compared to the 102-105 fields with SPIRE data alone. 
The numbers do not include the Witch Head Nebula, for which 
we have no PACS data and which therefore was excluded from 
this comparison. The background subtraction and the Sect. |A.2| 
corrections both decrease the mean value, but not significantly. 
Note that on physical grounds one could have expected the val¬ 
ues to increase with background subtraction (if dense material 
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has higher r(250yum)/r(/)) and to decrease with the inclusion 
of 160/mi channel (if the inclusion of shorter wavelengths in¬ 
creases estimated colour temperatures). The last column shows 
the difference relative to our default case. In this column we only 
list the common set of 81 fields where the error estimates are 
lower than 10% for all the five cases. On average, the changes 
in the r(250pm)/r(J) values are less than 1.0 units (lower than 
5%). The dispersion between different SPIRE analyses is lower 
than 2.0 units and somewhat higher when analyses of four and 
three bands are compared. 


Appendix B: Simulated NIR extinction maps 

In addition to photometric errors, the reliability of NIR optical 
depth estimates is mainly affected by the sampling provided by 
the background stars and the possible contamination by fore¬ 
ground stars and galaxies. 

We used 2MASS catalogue flags to eliminate most of the 
obvious galaxies. In addition to requiring a photometric quality 
corresponding to ph.qual in classes A-C, we excluded all point 
sources that were extended (flag extJcey is set) or were flagged 
with gaLcontam. These only remove part of the galaxies. The 
increased dispersion of intrinsic colours caused by galaxies is 
taken into account in the error estimates provided by the method 
NICER. Because our simulations use actual 2MASS data near 
the target fields, this extragalactic contamination is also auto¬ 
matically present in the simulations described below. 

The simulations are based on dust 250 pm optical depth 
maps derived from Herschel observations. Using only SPIRE 
data, we derive column densities at 25" resolution as a com¬ 
bination of 


r = r(500) + [r(350) - r(350 -> 500)], 


(B.l) 


where r(500) is calculated using 250, 350, and 500 pm maps 
at the lowest common resolution, r(350) is calculated similarly 
from 250/rni and 350/mi maps, and r(350 —> 500) is the lat- 
ter convolved to the re solution of the 500 pm observations (see 
Palmeiri m et al.||2013| ). Thus, the expression in square brack¬ 
ets describes structures that are seen at the resolution of 350/rni 
data (25"), but not at the resolution of 500 pm data (36"). The 
calculations assume a fixed value (3=2.0. The differences to the 
r(250//m) maps used in Sect. [3] are small and, furthermore, we 
only consider differences between these input maps and the val¬ 
ues recovered by NICER. On the other hand, we wish to retain 
the highest resolution possible (18" instead of 36") because the 
bias in NICER estimates is probably linked to the amount of 
small-scale structure. 

The Besancon model ( [Robin et al.|2003| ) was used to create 
a simulated catalogue of stars over a0.5° x0.5° area centred on 
each target field. The catalogue includes stellar distances and tl- 
band magnitudes, and together with the distance estimates listed 
in Table [I] this was converted into the probability that a star of 
given magnitude resides between the cloud and the observer. In 
the simulation, the corresponding fraction of stars was assumed 
to be located in front of the cloud. 

To simulate NIR observations, we used the same 2MASS 
data that were used to calculate the actual r(J) m aps of the fields. 
The stars in the reference area (see Table |E.1| ) were used to de¬ 
termine an empirical probability distribution of H-band magni¬ 
tudes and the dependence between the J, H, and Ks magnitudes 
and their uncertainties, as given in the 2MASS catalogue. This 
reference area may be affected by small amounts of absorption 
by diffuse dust, but gives a good approximation of the brightness 
distribution of stars that are unextincted by the main cloud. 


We simulated a uniform distribution of stars over the 
Herschel field, generating the magnitudes from the empirical 
H-band magnitude distribution and matching the average stellar 
density of the reference region. The J and Ks magnitudes of each 
star were generated using the J-H and H-Ks colours of a random 
star selected from the reference region. This ensures that the dis¬ 
tribution of intrinsic colours is realistic and that the simulations 
reproduce proper correlations (also in errors) between the bands. 

Based on the Besancon model, a fraction of stars was as¬ 
sumed to reside in front of the cloud and to be unaffected by ex¬ 
tinction. For the remaining stars, the line-of-sight optical depth 
in J band was calculated using the input column density map and 
a fixed ratio of r(250pm)/r(J) = 1.5 x 10 -3 . The optical depths 
in H and Ks bands then follow from the |Cardelli et al.l ( fl989| ) ex¬ 
tinction curve. The magnitudes were adjusted according to the 
line-of-sight optical depths. The magnitudes and their uncertain¬ 
ties in the reference area were used to calculate the typical pho¬ 
tometric uncertainties as a function of magnitude, and they were 
used in the NICER calculation. Because the intrinsic colours of 
the stars were generated based on observed stars, no additional 
scatter needed to be added for intrinsic colours. However, be¬ 
cause the extinction makes the stars fainter, the typical photo¬ 
metric errors increase as well. This was taken into account by 
adding normal distributed noise to the magnitudes, which corre¬ 
sponds to the difference in the typical uncertainties between the 
original and the extincted magnitudes. 

The simulated stellar catalogues were fed to the NICER al¬ 
gorithm to derive extinction maps with the same parameters as 
in Sect. |2.3| For each target field, one hundred realisations of the 
t(J) maps were calculated to obtain maps for the standard devi¬ 
ation and the bias of the estimated r(/) values. Figure [BT] shows 
one example. 


Appendix C: Simulated Herschel observations 

The estimation of the r(250 pm) bias is more uncertain than 
the estimation of r(J) bias. Because of line-of-sight tempera¬ 
ture variations, the colour temperature overestimates the mass- 
averaged dust temperature, which translates into too low esti¬ 
mates of r(250/im). The magnitude of the effect cannot be es¬ 
timated precisely because the line-of-sight temperature distribu¬ 
tion is unknown. Order of magnitude estimates can be obtained 
with radiative transfer modelling by making assumptions of the 
radiation field, dust properties, and the cloud structure. 

We carried out radiative transfer calculations individually 
for all the 116 fields. We assumed constant dust properties cor¬ 
responding to Milky Way dust with a selective extinction of 
Rv=5.5 ( |Draine|2003b|). The initial radiation field was assumed 
to correspond to the [Mathis et al.| ( |1983[ ) model of solar neigh¬ 
bourhood. The spectrum of the illuminating radiation has an im¬ 
pact on the temperature contrasts. We have no way to indepen¬ 
dently determine the shape of the spectrum of the radiation field. 
However, this typically remains a second-order effect, and the 
main effect, the level of the radiation field intensity, is part of the 
modelling. One exception are the possible stars that heat clouds 
from the inside and may locally have a strong effect. These are 
considered later in the analysis, but are not part of the simula¬ 
tions. 

For each field, we built a model that attempts to reproduce 
the 250-500 pm observations of that field. From the background- 
subtracted surface brightness maps we selected the central 30' x 
30' area. The model cloud had the same angular dimensions and 
was discretised onto a 181 3 cell grid. Each cell of the model 
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Fig.B.l. Field G300.86-9.00 as an example 
of r(/) bias maps estimated with simulated 
NICER observations. The frames show the op¬ 
tical depth derived from actual observations 
(frame a), average recovered extinction map in 
simulations (frame b), and the bias as the dif¬ 
ference between the output and input maps in 
the simulation (frame c). 


therefore corresponds to an angular size of 10" , but the lin¬ 
ear scale depends on the distance of the cloud. In the line-of- 
sight direction we assumed a Gaussian density distribution with 
a FWHM equal to 25% of the field size. The linear size again 
depends on the cloud distance because in more distant fields we 
are also probably concerned with larger structures. In the line- 
of-sight direction the density peak is always in the central plane 
of the model cube. This increases mutual shadowing of dense 
regions, which in reality can reside at different distances and 
increases the temperature contrasts in the model. On the other 
hand, for a field at 200 pc distance, the selected line-of-sight 
FWHM extent is ~0.4 pc, which is larger than the size of typ¬ 
ical cores. Therefore, it is difficult to say whether the selection 
of this particular line-of-sight density profile leads to an over- 
or underestimation of the final r(250yum) bias. These are usually 
second-order effects (Juvela et al.] ||2013| ) except for very dense 
clumps that can remain practically invisible in r(/) maps as well. 

We carried out radiative transfer calculations to produce syn¬ 
thetic surface brightness maps at Herschel wavelengths that were 
then convolved to the resolution of the observations. The ra¬ 
tio of observed and modelled 350 yum maps was used to adjust 
the column densities, applying the same multiplicative factor to 
all cells along the same line-of-sight. The intensity of the exter¬ 
nal radiation field was scaled based on the ratios between the 
observed and modelled 250/im and 500gm surface brightness. 
The aim is to also reproduce the average shape of the spectra. 
The full procedure was iterated until the model matched the 
observed 350 //m map at ~1% accuracy and the average ratio 
250 //m/500 /im were correct within the same tolerance. 

The final model takes into account the range of column den¬ 
sity, the morphology, and the radiation field intensity of a field. 
It is not necessarily a perfect match to all surface brightness 
maps, but is a good facsimile of the pixel-by-pixel column den¬ 
sity struct ure. We analysed the synthetic surface brightness maps 
as in Sect. |2.2.3| The comparison of the obtained r(250yum) es¬ 
timates and the true values known for the model cloud gives a 
30' x30' map of the r(250/rm) bias (resolution 40"). The remain¬ 
ing border areas usually have a low column density and therefore 
low bias. However, to extend bias estimates over the whole map 
area, we assigned values calculated using the average bias vs. 
column density relation estimated from central 30' x 30' area to 
the remaining pixels. 

Figure [CJ] shows one example, the surface brightness maps 
produced by the model and the bias in the r(250yum) values esti¬ 
mated from these synthetic observations. 


Appendix D: Additional checks of correlations 
between r(250yum) and r(J) 

In addition to the factors examined in Sect. |2.4[ we checked the 
importance of two additional factors, the technical implemen¬ 
tation of the least-squares fits, and the importance of internally 
heated regions. 


The total least-squares fits of Sect. |2.4| used the formal er¬ 
ror estimates of r(J) and r(250/im). The former were obtained 
from NICER routine, the latter were estimated with MCMC cal¬ 
culations starting with the assumption of 7% (SPIRE) or 15% 
(PACS) relative errors in the surface brightness data. If the corre¬ 
lation is poor, the estimated slope becomes sensitive to the error 
estimates. For example, if the true errors of r(250/im) were much 
larger, for example because of some artefacts in map making, 
the use of too low error estimates would increase the slope esti¬ 
mates. We checked this by repeating the analysis using twice the 
original r(250/im) error estimates. In an extreme case, we can 
ignore the error estimates altogether and perform unweighted 
least-squares fits. Based on the error estimates used, the true rel¬ 
ative uncertainty is significantly larger in r(J) than in r(250yum). 
Therefore, the unweighted least-squares fit probably underesti¬ 
mates the true slope. 

The third test concerns the internally heated regions in 
which because of strong temperature variations combined with 
compact, high column density clumps, both optical depth esti¬ 
mates are particularly uncertain. Furthermore, the estimates of 
r(250/rm) bias are probably incorrect in the same areas. This is 
caused by two factors. First, without the internal heating source, 
the models are unable to produce sufficient surface brightness 
values and result in very high column densities and thus high es¬ 
timates of the bias ( [Juvela et al.|2013| ). Second, in the real clump 
the internal heating may help to decrease the actual bias if the 
clump centre remains war m instead of being to o cold to be reg¬ 
istered in Hersc hel b ands ( [Malinen et al.|201 f[ ). We repeated the 
analysis of Sect. |2.4| by masking warm regions. We first masked 
all pixels for which the dust colour temperature was higher than 
20 K. The mask was then extended to cover areas in which after 
convolution to 180" resolution, the influence of the T > 20 K 
region was more than 10% of the convolved value. 

Figure |D.1| compares the result with the bias-corrected re¬ 
sults already shown in Fig. [4] It shows that the shape of the dis¬ 
tribution is not sensitive to the fitting procedure, nor is it signifi¬ 
cantly affected by the warm regions. 


Appendix E: Maps of r(250/im)/r(7) ratio 

In the least-squares fits of Sect. |3.4| the interesting parameters, 
k and C were independent of additive errors in the correlated 
quantities. The highest r(/) points were particularly important, 
both vi suall y and regarding the fitted parameters. As mentioned 
in Sect. |3.5| we also calculated maps of the r(250yum)/r(/) ratio. 
Because the column densities are typically low over most of the 
mapped area, the visual appearance of the maps is dominated by 
regions with low r(250yum) and r(/) values where, by definition, 
the results become sensitive to any zero-point mismatch. 

The maps were calculated by first correcting the r(250/rm) 
and r(/) maps for the bias that was estimated with modelling 
(see Sect. [B] and Sect. 0- The maps of r(250yum) were then 
convolved to the resolution of the r(J) map. A 2' wide region 
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Fig. C.l. Field G300.86-9.00 as an example of the r(250yum) bias estimation with radiative transfer modelling. The upper row 
shows the relative error between the model predicted surface brightness and the observations at 250 yum, 350/im, and 500 /am. The 
lower frames show the colour temperature and r(250/mi) maps calculated from the synthetic observations and the relative bias in 
r(250yum) obtained by comparison with the actual r(250yum) values in the model. 
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Fig. D.l. Comparison of Ar(250yum)/Ar(/) bias-corrected dis¬ 
tributions. In addition to the default case, derived distributions 
are shown for tests with larger r(250/im) error estimates, nor¬ 
mal unweighted least squares, and fits excluding data affected 
by regions with dust temperatures exceeding 20 K. 


near the Herschel map borders was masked because the con¬ 
volved values would be affected by data outside our map cover¬ 
age. Before calculating the ratio r(250yum)/r(/), we subtracted 
from both quantities the values in the reference regions that are 
listed in Table [T] A typical diameter of these reference areas is 
~ 6'. For r(250yum) the statistical error of the reference value is 
very small. The true uncertainty of the remaining r(250/mi) is 
completely dominated by systematic errors. Because the main 
purpose of the ratio maps is to determine variations in the 
r(250yum)/r(/) ratio, we are not very concerned with multiplica¬ 
tive errors. 

We expect the statistical errors to be more significant in r(/) 
and, because the variable is in the denominator, we need to mask 
areas with a low SN of r(/). NICER has provided error maps for 
r(/), but here we estimated the uncertainty using the following 
procedure: We took the data plotted in Fig. [8] selected 20% of 
the lowest r(/) points, and subtracted from them the prediction 
of the non-linear fit (red line in Fig. |8j. The uncertainty of the 
reference value, Ar(/), was calculated as the standard deviation 
of the residuals, scaled by the ratio of (90") 2 and the area of the 
reference region. This should be a very conservative estimate 
because it assumes that in Fig.[8]the scatter would be due to r(/) 
errors alone. 

The results are shown in Fig. |E.1| The first frames show the 
r(250//m) maps, the contours indicating the region with a SN 
of each parameter higher than one. The second frames show 
the calculated maps of r(250yum)/r(/). The regions where ei¬ 
ther parameter falls below SN=0.5 were masked. The remaining 
frames show the extreme cases corresponding to r(J) ± Ar(/) 
and r(250yum) ± 5t( 250yum), where cfr(250/mi) is the estimated 
error map (see Sect. |2. 2. 3] ). 
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Fig.E.l. Maps of r(250/im)/r(/) for the selected fields. The upper frames show r(250yum) (frame a) and the ratio r(250yum)/r(/). 
The lower frames show the lower (frame c) and upper (frame d) limits of r(250yum)/r(/) calculated as (r(250yum) + 
^r(250/im))/(r(/) - Ar(/)) and (r(250/rm) - ^r(250yum))/(r(/) + Ar(/)) where 6t( 250/rm) is the error map of r(250yum) and Ar(/) 
is the estimated uncertainty of the r(/) zero point. The areas in which the SN of either variable drops below 0.5 have been masked. 
In the first frame, the black contour corresponds to r(250yum) = cfr(250yum) and the dashed white contour to r(/) = Ar(/). 
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Table E.l. Coordinates and radii of the reference regions used to calculate the 
near-infrared optical depth maps. The last column lists the id numbers of the 
Herschel fields themselves. 


Field Reference area Radius Herschel Observation IDs 

RA(2000) DEC(2000) (') 


G0.02+18.02 

16 

43 

59.5 

-18 29 20.4 

16.8 

1342227737, 

1342216067, 

1342216066 

GO.49+11.38 

17 

01 

01.2 

-21 23 27.6 

10.1 

1342227645, 

1342216500, 

1342216499 

Gl.94+6.07 

17 

26 

46.8 

-23 17 09.6 

13.7 

1342216944, 

1342216588, 

1342216587 

G2.83+21.91 

16 

30 

59.3 

-13 27 46.8 

9.6 

1342227646, 

1342227040, 

1342227039 

G3.08+9.38 

17 

18 

49.9 

-20 40 01.2 

14.8 

1342252016, 

1342252015, 

1342251955 

G3.72+21.02 

16 

43 

08.9 

-13 19 37.2 

7.7 

1342227647, 

1342216502, 

1342216501 

G4.18+35.79 

15 

50 

46.1 

-05 04 37.2 

18.8 

1342215385, 

1342215384, 

1342213473 

G6.03+36.73 

15 

56 

12.2 

-02 02 45.6 

10.0 

1342204167, 

1342204166, 

1342203075 

G9.45+18.85 

17 

03 

07.7 

-10 09 54.0 

18.0 

1342227648, 

1342216504, 

1342216503 

G10.20+2.39 

17 

56 

04.1 

-18 09 07.2 

12.1 

1342217762, 

1342217761, 

1342216945 

G20.72+7.07 

18 

02 

13.7 

-07 34 51.6 

10.0 

1342216948, 

1342216592, 

1342216591 

G21.26+12.11 

17 

43 

37.9 

-03 47 27.6 

11.3 

1342216949, 

1342216590, 

1342216589 

G24.40+4.68 

18 

16 

02.4 

-04 43 19.2 

10.9 

1342217409, 

1342217408, 

1342216950 

G25.86+6.22 

18 

18 

31.7 

-03 17 56.4 

11.3 

1342217772, 

1342217771, 

1342216951 

G26.34+8.65 

18 

05 

08.9 

-01 01 22.8 

8.4 

1342216952, 

1342216594, 

1342216593 

G37.49+3.03 

18 

46 

05.3 

+05 53 42.0 

17.7 

1342219075, 

1342219074, 

1342216953 

G37.91+2.18 

18 

49 

47.3 

+06 09 46.8 

16.0 

1342229180, 

1342219077, 

1342219076 

G39.65+1.75 

18 

53 

57.8 

+07 30 46.8 

15.1 

1342219633, 

1342219079, 

1342219078 

G62.16-2.92 

20 

02 

27.4 

+23 30 36.0 

14.9 

1342219986, 

1342219043, 

1342219042 

G69.57-l.74 

20 

15 

55.7 

+30 40 19.2 

14.2 

1342220090, 

1342220089, 

1342219988 

G70.10-l.69 

20 

17 

18.0 

+31 06 46.8 

17.2 

1342220092, 

1342220091, 

1342219987 

G71.27-11.32 

20 

56 

03.6 

+27 36 25.2 

15.0 

1342244225, 

1342244224, 

1342244143 

G82.65-2.00 

20 

57 

19.2 

+40 46 30.0 

11.2 

1342222137, 

1342221287, 

1342221286, 1342219611 

G86.97-4.06 

21 

17 

01.4 

+42 51 00.0 

14.7 

1342222116, 

1342220280, 

1342220279, 1342219610 

G89.65-7.02 

21 

42 

49.9 

+44 06 32.4 

9.1 

1342220078, 

1342220077, 

1342213455 

G91.09-39.46 

23 

14 

05.0 

+ 17 15 18.0 

10.9 

1342221461, 

1342221285, 

1342221284 

G92.04+3.93 

20 

57 

44.9 

+53 18 28.8 

12.8 

1342220862, 

1342219031, 

1342219030 

G92.63-10.43 

22 

03 

29.5 

+41 39 25.2 

7.4 

1342220869, 

1342220080, 

1342220079 

G93.21+9.55 

20 

31 

49.7 

+57 36 39.6 

18.2 

1342220821, 

1342220820, 

1342213457 

G94.15+6.50 

20 

53 

42.5 

+56 08 31.2 

12.9 

1342220553, 

1342220552, 

1342213456 

G95.76+8.17 

20 

50 

34.3 

+59 01 55.2 

8.3 

1342244196, 

1342243763, 

1342243762 

G98.00+8.75 

21 

07 

10.8 

+59 13 26.4 

15.1 

1342220551, 

1342220550, 

1342220525 

G105.57+10.39 

21 

33 

44.9 

+67 18 28.8 

13.7 

1342219970, 

1342220819, 

1342220818 

G107.20+5.52 

22 

13 

59.8 

+63 38 09.6 

11.7 

1342187332, 

1342187331 


G108.28+16.68 

21 

07 

11.8 

+72 31 48.0 

13.5 

1342216072, 

1342216071, 

1342213452 

G109.18-37.59 

00 

07 

23.0 

+24 49 22.8 

9.0 

1342213511, 

1342213510, 

1342213494 

G109.80+2.70 

22 

49 

03.8 

+63 14 56.4 

16.5 

1342187655, 

1342187009, 

1342187008, 1342187007 

Gl 10.62-12.49 

23 

41 

20.9 

+47 50 16.8 

19.0 

1342222136, 

1342222110, 

1342222109 

Gl 10.89-2.78 

23 

17 

37.7 

+57 17 45.6 

13.0 

1342223343, 

1342223342, 

1342213453 

G110.80+14.16 

21 

49 

39.1 

+73 23 02.4 

12.4 

1342220622, 

1342220817, 

1342220816 

Glll.41-2.95 

23 

16 

59.8 

+57 02 52.8 

13.0 

1342223341, 

1342223340, 

1342213454 

Gl 15.93+9.47 

23 

20 

49.2 

+72 01 12.0 

8.8 

1342222598, 

1342220666, 

1342220665 

Gl 16.08-2.40 

00 

01 

23.8 

+58 57 28.8 

16.6 

1342222599, 

1342222158, 

1342222157 

G126.24-5.52 

01 

21 

49.9 

+57 04 19.2 

10.9 

1342224972, 

1342216074, 

1342216073 

G126.63+24.55 

04 

04 

11.0 

+84 55 55.2 

8.6 

1342243735, 

1342243734, 

1342219426, 1342219425, 







1342219424, 

1342199364 


G127.79+2.66 

01 

44 

44.9 

+65 20 38.4 

17.2 

1342216512, 

1342216511, 

1342203610 

G128.78-69.46 

00 

56 

41.0 

-06 16 12.0 

18.2 

1342246680, 

1342246679, 

1342246579 

G130.37+11.26 

02 

33 

16.1 

+71 50 20.4 

9.8 

1342218719, 

1342218718, 

1342216918 

G130.42-47.07 

01 

12 

09.8 

+ 14 48 43.2 

19.2 

1342213525, 

1342213524, 

1342213486 

G131.65+9.75 

02 

36 

24.2 

+71 39 21.6 

9.5 

1342243761, 

1342203609 

1342243760, 

1342223882, 1342218640, 

G132.12+8.95 

02 

43 

03.4 

+70 37 44.4 

11.4 

1342223881, 

1342223880, 

1342216919 

G139.60-3.06 

02 

47 

26.6 

+54 55 48.0 

16.5 

1342226625, 

1342216414, 

1342216413 

G141.25+34.37 

08 

37 

44.9 

+73 19 22.8 

13.3 

1342244197, 

1342243739, 

1342243738 

G149.67+3.56 

04 

15 

01.4 

+55 59 38.4 

14.1 

1342217533, 

1342217532, 

1342216920 

G150.47+3.93 

04 

22 

01.2 

+55 45 36.0 

17.1 

1342217531, 

1342217530, 

1342214702 

G151.45+3.95 

04 

31 

30.5 

+53 51 28.8 

9.9 

1342205043, 

1342205042, 

1342203607 
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Table E.l. continued. 


Field Reference area Radius Herschel Observations IDs 

RA(2000) DEC(2000) (') 


G154.08+5.23 

04 

52 

23.5 

+53 47 24.0 

18.1 

1342217529, 

1342217528, 

1342216921 

G155.80-14.24 

03 

40 

03.4 

+38 29 06.0 

13.0 

1342226626, 

1342224781, 

1342224780 

G157.08-8.68 

04 

02 

21.6 

+42 01 15.6 

13.1 

1342216416, 

1342216415, 

1342203615 

G157.92-2.28 

04 

28 

53.0 

+44 30 25.2 

13.7 

1342217527, 

1342217526, 

1342216923 

G159.23-34.51 

02 

59 

58.6 

+ 19 12 03.6 

13.8 

1342239264, 

1342239263 


G159.12-14.30 

03 

46 

43.0 

+34 50 09.6 

7.8 

1342226627, 

1342216419, 

1342216418 

G159.34+11.21 

05 

45 

49.0 

+51 29 31.2 

18.6 

1342218721, 

1342218720, 

1342216922 

G161.55-9.30 

04 

12 

05.5 

+37 11 20.4 

11.5 

1342205045, 

1342205044, 

1342203616 

G163.82-8.44 

04 

23 

32.6 

+36 22 48.0 

16.7 

1342205048, 

1342205047 


G164.71-5.64 

04 

43 

22.8 

+38 52 48.0 

14.5 

1342217525, 

1342217524, 

1342216925 

G167.20-8.69 

04 

34 

25.9 

+34 10 44.4 

6.6 

1342217392, 

1342217391, 

1342216926 

G168.85-10.19 

04 

33 

29.5 

+31 49 19.2 

14.1 

1342217523, 

1342217522, 

1342214701 

G171.35-38.28 

03 

19 

13.7 

+09 39 25.2 

12.8 

1342224970, 

1342224215, 

1342224214 

G173.43-5.44 

05 

07 

51.4 

+30 37 55.2 

17.3 

1342217507, 

1342217506, 

1342216927 

G174.22+2.58 

05 

42 

31.9 

+34 22 48.0 

9.9 

1342244179, 

1342243759, 

1342243758 

G176.27-2.09 

05 

24 

21.8 

+30 08 06.0 

10.3 

1342205197, 

1342205196, 

1342203617 

G181.84-18.46 

04 

47 

25.7 

+ 16 18 21.6 

10.3 

1342217469, 

1342217468, 

1342203625 

G188.24-12.97 

05 

14 

20.9 

+ 14 06 54.0 

11.9 

1342217455, 

1342217454, 

1342214700 

G189.51-10.41 

05 

28 

18.5 

+ 15 00 07.2 

10.9 

1342217457, 

1342217456, 

1342216930 

G195.74-2.29 

06 

13 

50.2 

+ 14 48 46.8 

17.9 

1342219409, 

1342219408, 

1342216931 

G198.58-9.10 

05 

55 

25.0 

+08 15 57.6 

17.9 

1342218702, 

1342218701, 

1342216932 

G202.23-3.38 

06 

18 

57.4 

+07 14 02.4 

10.2 

1342218773, 

1342218772, 

1342216933 

G202.02+2.85 

06 

42 

51.8 

+ 11 30 43.2 

12.2 

1342228372, 

1342228371, 

1342228342 

G203.42-8.29 

06 

07 

25.2 

+05 00 18.0 

12.6 

1342218777, 

1342218776, 

1342216935 

G205.06-6.04 

06 

19 

59.0 

+04 49 15.6 

7.6 

1342218775, 

1342218774, 

1342216934 

G206.33-25.94 

05 

04 

09.1 

-05 37 58.8 

9.3 

1342249237 



G210.90-36.55 

04 

31 

38.4 

-13 23 34.8 

9.7 

1342225213, 

1342225212, 

1342216940 

G212.07-15.21 

05 

59 

01.9 

-05 58 30.0 

14.0 

1342218783, 

1342218782, 

1342216936 

G215.37-3.04 

06 

44 

21.4 

-04 16 37.2 

15.7 

1342219957, 

1342219405, 

1342219404 

G215.44-16.38 

05 

59 

46.6 

-08 52 40.8 

19.9 

1342204306, 

1342204305, 

1342203631 

G216.76-2.58 

06 

45 

47.8 

-05 23 56.4 

12.0 

1342219956, 

1342219407, 

1342219406 

G218.06+2.12 

07 

11 

45.8 

-02 43 48.0 

8.5 

1342219958, 

1342220785, 

1342220784 

G219.36-9.71 

06 

28 

11.3 

-09 30 07.2 

15.6 

1342227708, 

1342219403, 

1342219402 

G219.29-9.25 

06 

28 

20.9 

-09 27 43.2 

16.5 

1342227707, 

1342219401, 

1342219400 

G227.95-2.98 

07 

04 

29.5 

-15 27 39.6 

10.9 

1342219982, 

1342220901, 

1342220900 

G247.55-12.27 

07 

13 

20.6 

-36 48 00.0 

12.7 

1342222831, 

1342220781, 

1342220780 

G253.71+1.93 

08 

25 

07.0 

_33 44 45.6 

10.8 

1342222828, 

1342220783, 

1342220782 

G255.33-4.88 

07 

58 

19.9 

-39 40 51.6 

16.7 

1342222830, 

1342220779, 

1342220778 

G258.90-4.10 

08 

10 

17.0 

-42 19 15.6 

10.6 

1342222896, 

1342220777, 

1342220776 

G265.04+6.08 

09 

21 

18.7 

-40 14 38.4 

7.8 

1342222845, 

1342220309, 

1342220308 

G265.60-5.82 

08 

32 

48.5 

-47 47 06.0 

12.0 

1342222829, 

1342220775, 

1342220774 

G268.21+2.02 

09 

17 

09.8 

-44 53 49.2 

17.7 

1342222827, 

1342221275, 

1342221274 

G271.06+4.84 

09 

41 

12.0 

-45 28 51.6 

9.9 

1342222895, 

1342221270, 

1342221269 

G271.51+5.14 

09 

44 

13.9 

-45 54 39.6 

10.2 

1342222894, 

1342221268, 

1342221267 

G276.78+1.75 

09 

54 

27.1 

-51 06 21.6 

9.8 

1342198593, 

1342198592 


G298.31-13.05 

11 

26 

36.0 

-75 02 52.8 

11.7 

1342223601, 

1342223600, 

1342216941 

G299.57+5.61 

12 

22 

23.3 

-56 28 33.6 

23.5 

1342223605, 

1342223604, 

1342223259 

G300.61-3.13 

12 

27 

23.0 

-66 32 20.4 

15.1 

1342223603, 

1342223602, 

1342213480 

G300.86-9.00 

12 

14 

18.7 

-72 08 02.4 

11.4 

1342188162, 

1342188101, 

1342188100, 1342188099 

G315.88-21.44 

17 

06 

03.4 

-77 16 37.2 

14.6 

1342218738, 

1342218737, 

1342216942 

G320.84+5.09 

14 

55 

25.2 

-52 38 56.4 

13.2 

1342227725, 

1342215598, 

1342215597 

G325.54+5.82 

15 

16 

29.0 

-49 30 54.0 

9.8 

1342227693, 

1342225000, 

1342224999 

G332.70+6.77 

15 

47 

46.1 

-44 43 04.8 

15.5 

1342227692, 

1342226700, 

1342226699 

G334.65+2.67 

16 

10 

07.0 

-46 25 26.4 

13.5 

1342216490, 

1342216489, 

1342214758 

G339.22-6.02 

17 

15 

55.9 

-50 25 01.2 

9.3 

1342216582, 

1342216581, 

1342214757 

G341.18+6.51 

16 

21 

30.2 

-39 15 07.2 

9.4 

1342227691, 

1342216492, 

1342216491 

G343.64-2.31 

17 

14 

32.4 

-44 28 40.8 

17.3 

1342216943, 

1342216584, 

1342216583 

G344.77+7.58 

16 

29 

49.0 

-37 18 21.6 

13.7 

1342227690, 

1342216496, 

1342216495 

G345.39-3.97 

17 

26 

34.1 

-43 44 13.2 

20.6 

1342252027, 

1342251963 


G358.96+36.75 

15 

36 

34.1 

-07 12 54.0 

12.2 

1342204169, 

1342204168, 

1342202205 
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